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Tensor networks, and in particular Projected Entangled Pair States (PEPS), are a powerful tool 
for the study of quantum many body physics, thanks to both their built-in ability of classifying and 
studying symmetries, and the efficient numerical calculations they allow. In this work, we introduce 
a way to extend the set of symmetric PEPS in order to include local gauge invariance and investigate 
lattice gauge theories with fermionic matter. To this purpose, we provide as a case study and first 
example, the construction of a fermionic PEPS, based on Gaussian schemes, invariant under both 
global and local U{1) gauge transformations. The obtained states correspond to a truncated U{1) 
lattice gauge theory in 2 + 1 dimensions, involving both the gauge field and fermionic matter. For 
the global symmetry (pure fermionic) case, these PEPS can be studied in terms of spinless fermions 
subject to a p-wave superconducting pairing. For the local symmetry (fermions and gauge fields) 
case, we find confined and deconfined phases in the pure gauge limit, and we discuss the screening 
properties of the phases arising in the presence of dynamical matter. 
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I. INTRODUCTION 

Within the framework of the standard model of particle physics, the three fundamental forces are described by 
gauge bosons, which are the excitations of gauge fields. Gauge fields are vector fields, which manifest a very special 
local continuous symmetry, called local gauge invariance. This symmetry gives the matter fields gauge charges, and 
its local nature induces local interactions of the gauge currents with the charged matter. The conservation of local 
charges, manifested by local constraints which are extensions of the well-known Gauss law from electrodynamics, 
implies a very rich, complicated structure of the Hilbert space of quantum gauge theories, dictated by superselection 
rules governed by these local charges. This makes such theories, in general, very challenging and difficult to solve - 
just as much as they are interesting and important for the description of nature. 

Described within the framework of quantum field theory, gauge theories come along with a very important com¬ 
putational tool: perturbation theory and its Feynman diagrams. However, despite the great success and accuracy 
achieved for Quantum Electrodynamics (QED) with perturbative methods, they apply only partially to QCD. Unlike 
QED, which is the Abelian gauge theory associated with the group 17(1), QCD is a non-Abelian, SU{3) gauge theory, 
which makes it behave in a completely different manner: due to an important property of such non-Abelian theories, 
Asymptotic Freedom 0], the strong coupling constant flows to zero for high energies (or short distances) - allowing, 
therefore, for perturbative calculations in these scales, such as within the nuclei (e.g., the parton model, or Bjorken’s 
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scaling mm)- On the other hand, at low energies or large distances, the coupling constant is strong and perturbative 
physics is impossible; this may be seen as both the cause and the effect of Quark Confinement 015], the phenomenon 
responsible for holding quarks bound together into hadrons, and for the absence of free quarks in the spectrum of the 
theory. 

This has significant implications on the study of the theory, and, indeed, over the years many non-pertubative 
techniques have been developed and applied for the study of QCD, and non-Abelian gauge theories in general. One 
of them, perhaps the most fruitful, is lattice gauge theory (LGT) 0|MH]) in which either spacetime, or space, is 
discretized, allowing either for a regularization of the theory for analytical purposes, or very efficient and fruitful 
numerical (Monte Carlo) calculations. While having a great success with many different types of calculations and 
predictions (e.g., low-energy hadronic spectrum |^, among others), Monte Carlo calculations are problematic in 
some cases: first, with fermions with a finite chemical potential (required, for example, for the phases of color 
superconductivity and quark-gluon plasma [IQIIII]), due to the computationally hard sign problem [12) . and second, 
as the calculations are carried out in Euclidean spacetime, real-time dynamics in Minkowski spacetime cannot be 
achieved (see, on the other hand, the recent works [Hill]). 

A complementary way of overcoming the computational difficulties may be the use of tensor network techniques, 
or tensor network states (TNs), and in particular Matrix Product States (MPS) [TS] and Projected Entangled Pair 
States (PEPS) [TBUIS] . One may, for example, use TN variational techniques, in which a TN state with variational 
parameters is used as an ansatz for the ground state of a given Hamiltonian, as well as calculate dynamics of such 
states in very efficient ways, exploiting methods like DMRC (Density Matrix Renormalization Croup) [12]. This 
approach has been recently applied with MPS for 1-1-1 dimensional lattice gauge theories, either Abelian or non- 
Abelian, and used for the study of their spectrum, dynamics (including string-breaking) and finite temperature effects 
PD1 - E2] . Furthermore also tensor renormalization group techniques have been recently applied to the study of such 
models [50] • The MPS studies have shown many of the static and dynamic properties of some well known theories 
(such as the Schwinger model |3T1|32|, for example), and allowed to reach better precision than analogous Monte Carlo 
calculations and to perform dynamical simulations, holding the promise for even better accuracy and computational 
possibilities. 

However, the great computational power is not the only reason for TNs to be candidates for the study of gauge 
theories. In a somewhat change of paradigm, one may describe a physical system from the point of view of its most 
representative states, instead of starting from its Lagrangian or Hamiltonian formulation. Tensor networks are, indeed, 
well suited to define families of states, as functions of a set of variational parameters, which fulfill precise symmetry 
constraints. Therefore they provide a natural way to encode all the symmetries of a system [33] and to describe its 
possible thermodynamical phases in terms of representative states, allowing to investigate the main physical properties 
within the universality class of the problem under scrutiny |34j . Starting from the tensor network construction, it is 
also possible to show that such states constitute the ground states of local parent Hamiltonians. These Hamiltonians 
may be explicitly derived in the simplest cases, and offer, as a function of the variational parameters, suitable examples 
to study the properties of the thermodynamical phases in a certain universality class. Previous works in this direction 
include two-dimensional PEPS schemes for lattice pure-gauge theories (without dynamical matter) |35j . as well as a 
general framework useful to the study of lattice gauge theories with bosonic matter |36j . 

The next reasonable step is to consider lattice gauge theories with fermionic dynamical matter, as in the case of 
high energy physics theories, for example. This paper addresses precisely this problem; more specifically, we analyze 
how one could utilize PEPS for the study of lattice gauge theories with dynamical matter. Could one classify locally 
gauge invariant states using PEPS, in a way that allows, eventually, to study the theories described by the Hilbert 
spaces to which they belong? 

For that purpose, we systematically construct, in this paper, fermionic PEPS (fPEPS) [37] in 2-1-1 dimensions, which 
have both local gauge symmetry and fundamental physical symmetries - rotation, translation and charge conjugation. 
To demonstrate the strength of fPEPS for studying such theories, we consider, as a case study, a truncated compact 
[/(I) gauge theory. Although simple, such states encode all the crucial ingredients for our demonstration: fermionic 
matter with bosonic gauge fields; nontrivial manifestation of the spatial symmetries, due to the fermions; and a rich, 
interesting phase diagram. 

We shall hereby show that using PEPS, one may capture the symmetry properties of a gauge theory, which are 
essentially the ones which define it, as can be deduced from its name; that both global and local symmetries may be 
manifested by PEPS - i.e., that one may use this method to treat both matter and gauge fields; and that PEPS allow 
to study the phase diagram of a gauge theory, and in some cases, using standard techniques, also to derive parent 
Hamiltonians. We shall emphasize, on the other hand, that the goal is not to study a compact 17(1) lattice gauge 
theory in 2 -|- 1 dimensions, but rather to show that PEPS may be used for the study of gauge theories, once the 
formalism we introduce is combined with efficient numerical methods. And, eventually, since PEPS allow us to find 
local parent Hamiltonians, most likely within the universality class of the model in question |34j . one may deduce that 
even if the parent Hamiltonian of a state in question, using the methods presented below, is not the one of the desired 
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lattice gauge theory, it is highly probable that this parent Hamiltonian will be in the same universality class, which 
means, that the two Hamiltonians would share many features. As the states in study are exact, i.e. exact ground 
statee of parent Hamiltonians, and one can perfoem numerical calculations, this could be used as a ”lab” for other 
theoretical methods used in high energy physics. 

Another possible avenue of exploring lattice gauge theories, which suggests a way of overcoming these difficulties, 
is quantum simulation. In recent years, many proposals have been made, for the mapping of lattice gauge theories. 
Abelian and non-Abelian, to atomic and optical systems (such as ultracold atoms in optical lattices, for example) 
[381139| . These systems, called quantum simulators, may be built in the laboratory and serve as quantum computers 
especially tailored for the purpose of lattice gauge theory calculations. Due to experimental requirements, one mostly 
has to approximate the simulated model by another one, with truncated local Hilbert spaces for the gauge degrees 
of freedom [351 l38H43j . A very important issue is the evaluation of the truncated approximation, which may also be 
done by the use of tensor network states [53J [55] . 


II. OUTLINE 


The work presented in this paper is organized as follows. First, in section m we introduce a class of states for 
staggered fermions |44| on a two dimensional spatial lattice, constructed as fermionic Gaussian PEPS |37|. These 
states will depend on a set of three parameters we shall introduce, in a way that guarantees both the spatial symmetries 
of translation and rotation invariance, and a global U{1) symmetry. As PEPS, these states are the ground states of 
local parent Hamiltonians. By constructing them as Gaussian PEPS, we will be able to easily derive these quadratic 
Hamiltonians, which will be, in our case, BdG Hamiltonians of p-wave superconductors. We will study the phase 
diagram of these states as a function of their parameters, and see that they exhibit gapped phases with a strong 
p-wave pairing, separated by gapless lines with either strong or weak pairing. 

Then, in section jrv] we will introduce new degrees of freedom, corresponding to a lattice Abelian gauge field, as in 
compact QED jd] |7| , but with finite (truncated) local Hilbert spaces |45l447j . This will allow us to gauge the U{1) 
symmetry of the states and make it local. Although these states, as these of an interacting theory, will no longer 
be Gaussian, they can still be parameterized exactly as the fermionic states with the global symmetry, manifesting 
the rotational invariance again and extending the translation invariance to a charge conjugation symmetry. From the 
general theory of PEPS [34] we know that, also in this case, the state [■i/'b) can be described as the ground state of 
a local parent Hamiltonian. There is a standard method for constructing the parent Hamiltonian, but deriving its 
explicit form is a lengthy procedure, out of the scope of our work which aims at describing the states. 

After the construction of the set of locally gauge-invariant states, we will be able to determine a phase diagram, 
as a function of the parameters of the states, by monitoring the transfer matrix of the PEPS as shall be explained. 
Furthermore, the tensor network approach allows us to calculate the expectation value of several important observables 
for these states, as, for example, the Wilson loops. We will find that the obtained states exhibit a rich variety of 
different behaviors, consistently with the expectations for the thermodynamical phases of a gauge theory - such as 
phases which confine and do not confine static charges in the case of pure-gauge states, and phases with different 
screening properties for states which involve, besides the gauge field, dynamical fermions. 

We do not assume familiarity of the reader with PEPS or tensor network states, therefore the paper is written and 
structured such that the results and their significance are clear even for people who are not experts in the field. The 
appendices include detailed derivations and proofs in PEPS language, supporting the main paper, and are written 
in a self-contained way, aimed both at the expert and the non-expert readers. Throughout the paper, the Einstein 
summation convention (summation on double indices) is assumed. 


III. GAUSSIAN FPEPS: GLOBAL SYMMETRY 


A. PEPS construction of globally invariant Gaussian states 


Hereby we shall describe the construction of globally invariant Gaussian states which shall fulfill several symmetries, 
yet to be classified. We assume no acquaintance of the reader with PEPS, and will thus describe in detail the process 
of constructing the state. 
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FIG. 1. The Hilbert space on a vertex (composing the fiducial state); a single physical fermion -i/), with eight virtual fermions 
surrounding it, two on each of the edges intersecting at the vertex. 


1. The PEPS construction 

Consider a two dimensional lattice (corresponding to a 2 + 1 dimensional spacetime), whose vertices (sites) are 
denoted by x G Z^, with unit vectors ei^ 2 - On each vertex x one defines a fermionic Fock space TLx [48], with a single 
physical mode annihilated by ipx, such that its vacuum state jUp (x)) satisfies 


tpx (x)) = 0. (1) 

We decompose this lattice into two sublattices. The even one (were the indices Xi and X 2 have the same parity) may 
be occupied by particles, and the odd by antiparticles. We define their charges accordingly, 

Qx — ; ( 2 ) 

where Sx = (—1)“^^®^. This is a definition of a staggered charge for staggered fermions. A particle-hole transformation, 
which we shall not carry out, maps the model into the Fock space of the Kogut-Susskind staggered fermions [44] . 
The staggering procedure allows us to analyze the system as a discretization of a theory with two-component spinors, 
staggered between the two sublattices. Such theory can be considered as embedded in the 3-1-1 dimensional theory 
by Susskind, for example, by taking the restriction X 3 = 0. In such a restriction only two spinorial components are 
required, which are, in the convention we use, the first and the fourth in the 3-1-1 dimensional construction of |44j 
(see also [1511301 for the Hamiltonian formulation in 2 -f 1 dimensions). 

Before considering the state of the entire lattice, we shall focus on the state of a single vertex x, called the fiducial 
state \F (x)), involving the already introduced physical mode '0^^, and other modes, called the virtual ones, which shall 
now be introduced and later will be traced out, while tailoring the vertices to each other, constructing the desired 
state \'ijj). 

On each vertex we define eight virtual modes, located on the edges of the links intersecting at the vertex. They are 
all fermionic, denoted by the annihilation operators l±,r±,u±,d±, which stand for “left”, “right”, “up” and “down” 
(see figure Their vacuum state is denoted by and the vacuum state of all the local fermions (both physical 
and virtual) by |H). 

If one concatenates all the physical and virtual creation operators on a given vertex to a vector of operators a|, the 
most general fiducial state, which is Gaussian, takes the form m 

|F(x))=A(x)|H(x)) (3) 


with 


A = exp 




(4) 


Here and in the following, while working at a given vertex, we will neglect the position index x. 

The next step is to connect the fiducial states lying on the lattice vertices among themselves, and project out the 
virtual degrees of freedom, for the creation of a physical state \if) for the entire lattice. To this purpose, we project 
the virtual states on both sides of a bond into joint entangled states (see figure]^, which allow to connect the whole 
network of fiducial states, hence the name Projected Entangled Pair States (PEPS). In particular, for each horizontal 
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FIG. 2. The ellipses denote the entangled state \H), onto which the virtual fermions are projected, for the contraction of the 
fiducial states into the PEPS. 


and vertical bond, we project both the positive and negative virtual modes on the following entangled bond states: 

^ exp (r; ^ ) exp (’^-.x ) 

(5) 


|i?x) = 2 exp exp 

l^^x) = 2 exp (w^+.xrf+,x+ej exp (uL.xt^l,x+e,) I^T,x) 

where \i^H,v) are the respective vacua on the links. From these states, one constructs the projection operators 

uj (x) = |iJx> {H^\ , ?7 (x) = |t^x) (f^xl • (6) 

Finally, the physical PEPS is given by: 

IV'(r)) = J|w(x)? 7 (x) A(x) |II). (7) 

X 

where |f 2 ) is the global vacuum, while jn^,) is the one of the virtual fermions only. 

2. The required symmetries 

Our goal is to construct a physical state \'ip {{ti})) of the physical modes, parameterized by a set of variational 
parameters {ti}, which satisfies the following symmetries: 

1. Translational invariance. Define the unitary operators Ut, such that 

Ut (ei) (ei) = 

Ut ( 62 ) V'iC/r ( 62 ) = 

The state {ip ({^i})) is translationally invariant if and only if 

Ut (en) IpJ {{ti})) = \ip {{U})) , n = l,2 


( 8 ) 


(9) 

(in general, invariance might also involve a change of the state by a global phase, but we shall neglect this 
option here and in the following symmetries; this makes sense in this case, since one naturally uses PEPS for 
zero momentum states; otherwise, different tensors are required, e.g. m)- 

2 . Rotational invariance. As we are dealing with a model defined on a square lattice, only C 4 rotations are relevant. 
We denote a counter-clockwise 7 r /2 rotation by A, i.e. 


X = {Xi,X2) —S> Ax = (-X2, Xi) . 

Denote by Up the quantum unitary operator which rotates the physical modes: 

Up'tp:x.Up = "^pV^Tlx 

where, in general, \r]p\ = 1. We will make, however, the choice of 


( 10 ) 


( 11 ) 


Vp = e*"/" 


(12) 
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which satisfies the set of physical requirements imposed by the staggered fermion discretization |44j of the Dirac 
theory [52]. A rotationally-invariant state has to satisfy 

Up{A)\^{{U})) = \i^{{U})) (13) 

(again, we are not allowing for a global phase). 

3. Global U{1) Invariance. We wish our state to be invariant under a global U{1) gauge transformation. Due to 
the staggering, this is the transformation generated by 

Go = (14) 

X X 

i.e., under a gauge transformation with the parameter (p, 

ipl -)■ (15) 

A globally gauge invariant state satisfies: 

e*<^^«|V^({tJ)) = |^({tJ)). (16) 

While these symmetries can be met by states which are either Gaussian or not, we shall use Gaussian PEPS as 
described above which shall satisfy these symmetries by construction. We will see how to obtain these symmetries 
first on the virtual level of the fiducial states, and then, using the PEPS projectors, to make them a real physical 
symmetry. 


3. Imposing global invariance 

On each of the tensor network links we define a virtual electric fields 

Ea = Q!^(q!+ — a^_a- 

(where a = l,r,u,d), with eigenvalues —1,0,1. We further define a local virtual Gauss operator, 

Go = divA — Q = Er + Eu — Ei — Ea — Q = 

r\_r+ — rir_ + u\_u+ — u^_u- — ll_l+ + l^_l- — d^^d^ + d^_d- — Sy^ip^ip 

We wish to define a fiducial state of the physical and virtual fermions, \F). This will be an eigenstate of the Gauss 
operator. 


(17) 

(18) 


Go\F)=q\F) 


(19) 


- a virtual Gauss law, with a static charge q. As we shall see, this will ensure, once we propely construct the PEPS, 
that we will eventually have a globally invariant state. 

We can decompose the physical and virtual modes into two sets with respect to this Gauss law: modes whose 


number operators appear in (18) with a negative sign will be called negative modes, and modes with a positive sign 


- positive modes. We shall denote the virtual negative modes by {ai}j^_^ = {l+,r-,U-,d+}, and the virtual positive 
modes by {bi}^^i = {/_,r_|_,rt+,d_}. The physical mode is negative on the even sublattice and positive on the odd 
one, and shall be denoted, in this notation, by either oq or bo for even or odd vertices, respectively. 

Note that, in general, one could introduce further virtual fermions, and update the Gauss law accordingly. In the 
general case, one may consider iV„ negative modes and Np positive ones; we will focus on the case \Nn — Ap| = 1, 
with two fermions per bond. 

We wish to find the most general Gaussian state \F) which is an eigenstate of Go, and of the form ([^, with two 
fermions per bond. Gaussian states are fully characterized by their covariance matrix T, and thus one may exploit 
the covariance matrix for their parametrization, as done in [^ Below we shall describe an equivalent way for their 
parametrization, which we shall utilize next. 



Statement 1. ; The most general gauge invariant Gaussian fiducial state |i^) with no static charges takes the form 


|F (x)) = A (x) (x)) (x)) = A (x) |fJ (x)) 


( 20 ) 


using the Gaussian operator 


exp I 1 > X 


A(x)= 




exp I 1 ; ^ 


( 21 ) 


i = For even vertices, Nn = 5, Np = 4, and the opposite for odd vertices. Thus, in general, 

before demanding any other symmetries, we obtain that the state depends on 20 complex parameters - the elements 
of the T matrix. 


Proof. The most general Gaussian fiducial state (even or odd) is constructed out of the Gaussian operator [33 


A = exp j 


ya|a] 


( 22 ) 


with i,j = 1,..., 9. The fiducial state is an eigenstate of Gq, if and only if 

Let us consider an infinitesimal transformation, with 0 = e <C 1. The transformation reads 


A^ Jfi 


exp 


E4 


(«! 


)a] + ie 


Go, ala] 


(23) 


(24) 


The desired result is obtained if the commutators 


Go,a|a] 


are c-numbers for any i,j. Go consists of number 


operators with different sign. If a| = a| and a] = 6], the commutator vanishes. If, instead, both the creation operators 
are either positive or negative, the commutator results in a term which is quadratic in the creation operators, and 
these terms cannot cancel each other. Thus we conclude that Tij is nonzero only if a| and a] create opposite charges, 
therefore a| = a|, a] = 5] (or the other way around). We also obtain that <; = 0. □ 

One can easily verify that the state Q is gauge-invariant: the global transformation we wish to apply on the state 


IS 


U = exp ; 


and thanks to statement [3 we know that 

\F (x)) = |F (x)) = |F (x)) = Uf, (x) Ulf (x) u'^ (x) Up (x) |F (x)) 

(where Ufi = ), and 

IJU yit ijdt ^ jjd ^ 

On the other hand, the projectors are invariant under these virtual operations, since 

l^x) = |i?x) , l^x) = |4"x) 

and thus one obtains that the state is, indeed, globally invariant: 


(25) 


(26) 


(27) 


(28) 


U\i;{T)) = \f;{T)). 


(29) 
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4- Virtual phase symmetries 


The construction of the PEPS is not unique, and includes a virtual symmetry, which is due to a redundancy in the 
definition of the virtual modes of such tensor networks m- In our construction, we have virtual symmetries which 
are defined by the following phase transformations on the virtual modes: 


Us {a, P, 1,5) a] Ul {a, /?, 7 , <5) = SAij {a, /3, 7 ,6 ) a] 


Us (a, 13, 1 ,5)h\ul (a, /3, 7 , 5) = Sbzj {a, /3, 7 ,6 ) h] 


(30) 


where 




for an even vertex, and 


^“4 = 


/ 1 0 0 0 0 \ 

' 0 e*“ 0 0 o' 

0 0 0 0 

0 0 0 e*T' 0 

V 0 0 0 0 6*“^ y 


e“ 0 0 0 \ 

0 0 0 
0 0 0 
Vo 0 0 y 


S% = 


e ' 

0 

0 

V 0 


0 

„ — 20 ! 


0 

0 


0 \ 

0 


0 e 

0 0 e 


0 

— 27 


(31) 


CO _ 
— 


1 0 
0 

0 0 
0 0 
Vo 0 


0 

0 

0 

„—iS 


0 0 e 

for an odd one. Note that the projectors are invariant under this transformation: 

Us (a, P, 1,3) ujUI (a, /3, 7 ,5) = w 
Us (a, P, 1, <5) rtUl {a, /?, -f,6) =r) 


^ \ 

0 

0 

0 

— 27 


(32) 


(33) 


The application of this transformation on the virtual modes does not affect the physical state. Let us define \tp{T)) 
as the state created using the operators A = exp ^(and similarly for odd vertices) where 


r = Sf {a,P,-f,S)TS%{a,p,^,d) , 
f‘^ = Sf {a,P,^,S)TS°A{a,P,-f,S) 

for an arbitrary choice of the phases a, P, 7 ,6 . Since T®’° result from A (x) = Us (x) A (x) Ug (x), we obtain 
|^/>(T)) = (rj„| ]V[w (x) 7 ] (x) A (x) |f 2 „) IVtp) = ( 12^1 ]V[a; (x) 7 (x) Us (x) A (x) Ul (x) |f 2 „) \ilp) = 

X X 

(rj„| ]V[C /5 (x) UJ (x) T] (x) Us (x) A (x) |f2„) irip) = \tp (T)) 


(34) 


(35) 


where we used the invariance of the virtual vacuum under the transformation. Therefore, we conclude that the 
transformations S are just virtual symmetries of the state |'0(r)), or consequences of a redundant description, which 
will allow us to reduce, eventually, the number of relevant parameters (especially the phases) of the parametrization 
matrix T. 


5. Imposing the Physical Symmetries 


So far. I'll) (T)) fulfills the global symmetry requirement, and depends on 20 complex parameters - the matrix 
elements of T. Demanding the other two required physical symmetries discussed in the beginning of this section, and 
using the virtual symmetry presented above, we shall now reduce the number of degrees of freedom further. 

We begin with rotational invariance. It follows from the definition of A (x ) in the form of a fermionic Gaussian op¬ 
erator, that such invariance cannot induce a global phase, consistently with ( |13[ ). We now proceed to the construction 
of the rotation transformation presented in ( 10 ),( 111 ,( 12 ). 
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FIG. 3. Rotation of a single fiducial state. 


In order to parameterize T such that rotational invariance is fulfilled, we have to define a rotation transformation 
for the virtual modes, in terms of a quantum unitary operator Un (/I), as follows: 


4 (x) — (yl) (x) Ur (yl) = r]uvl^ 

4 (x) —^ (yl) 4 (x) Ur (A) = r]ddi (yl"^x) 
(x) —^ u}^ (yl) vJ^ (x) Ur {A) = (vl"^x) 

dt (x) —C/fl (^) (x) Ur (A) = 77/Z| (yl"^x) 


(36) 


where r^/, Vd are phases we will specify in the following, and we described a clockwise rotation for convenience. 

This rotation preserves the Gauss law orientation (see Fig. [^. In particular, we impose rjuVd = I and 77 ^. 77 / = —I, in 
order to obtain the expected rotation transformation for the projectors: 


(37) 


Ur (yl) w (x) C/jj (yl) = r] (ylx) 

Ur (yl) 77 (x) i/Jj (yl) = w (ylx) 

Finally, to fulfill the rotational invariance, we adopt the parametrization of T specified by the following statement: 
Statement 2. The PEPS \ijj (T)) is invariant under the rotation transformation Up, if 

T = —y X —w z (38) 


t 

-Vp^t 

-Vp^t 

-% 

X 

y 

z 

w 

-y 

X 

—w 

z 

—w 

z 

X 

y 

—z 

—w 

-y 

X 


- reducing the number of free parameters to five complex ones in total. 

Proof: If we parametrize T (and therefore A) such that the entire fiducial state is rotated as a whole piece, with 
no change to its internal structure. 


(39) 


Up (A) Ur (vl) yl (x) u; (A) ul {A) = yl (^x) , 
then we get the required symmetry by using Eq. ( [T7| : 

Up (yl) IV') = {^v\ (x) 77 (x) Up (yl) A (x) ul (yl) |ll„) Iflp) 

X 

= {^v\ (x) 77 (x) c/jj (yl) A (ylx) Ur (yl) |n„) \cip) 

X 

= n Ur ^ w d W (A) yl (^x) \n„) |Gp) 

X 

= {^v\ (Tlx) uj (ylx) A (ylx) \ny) irip) = |V'). 


Therefore, the only remaining task is to parameterize T accordingly. The transformation Up (yl) Ur (yl) acts on the 
creation operators as 


(40) 


Ul (A) Ul (vl) oIUr (^) Up (^) = RA,ijal 
Ul (T) uj, (^) blUR (A) Up (A) = 


(41) 
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with the matrices 


/ 

Vp 

0 

0 

0 

0 

\ 


0 

0 

0 

rju 

0 


0 

0 

0 

0 

Vd 



0 

0 

Vr 

0 

0 


\ 

0 

Vi 

0 

0 

0 

/ 


(recall that rjd = Vu Vi = ~Vr) 


Rb 


0 0 ? 7 „ 0 \ 

0 0 0 r]d \ 

0 r]r 0 0 I 

m 0 0 0 / 


(42) 


(43) 


for an even vertex; for odd vertices, the two matrices have to be exchanged. Thus, the rotational symmetry is 
equivalent to imposing 


T=R\TRb 


(44) 


(where j labels the transposed matrix) for even vertices, and an analogous condition for odd vertices is obtained upon 
exchanging A and B. This is satisfied if and only if ry)) = — 1, consistently with our physical demand rjp = . 

The most general parametrization satisfying (44) takes the form: 


t 

-Tlr'nuVp'^t 

-VuVp 

-VrVp \ 

X 

y 

z 

w 

-y 

vlriix 

-vlw 

VrZ 

-rj-^rjuW 

VrVuZ 


y 

-Vu^VrZ 

-VrVuW 

-y 

■qlx / 


(45) 


Some phases of T^, though, can be eliminated using the fact that the state is left invariant under the transformation 
(30). By properly choosing the eight phases a,/3, 7 , i5, d,/3, 7 , <5. One may find out that the parametrization (45) 
is left invariant only if some requirements are fulfilled by the eight phases a,/3, 7 ,d,/3, 7 , <5, namely, one is only 
left with a single phase, out of which all the others may be expressed. Then, using the transformation matrices 
S% = Sq = 1 (B S, = S% = S, with 


/ 0 0 0 \ 



we find the equivalent parametrization 


T = 



-VpWVrVut 

-Vp" y'llr'nut 

-Vp 

VrVuX 

y 

VrZ 

VuW 

-y 

rjrtluX 

-VuW 

VrZ 

-riuW 

VrZ 

VrVuX 

y 

-VrZ 

-VuW 

-y 

VrVuX 


\ 

/ 


(46) 


(47) 


for both even and odd vertices. From this form, it is clearly seen that one may absorb the phases ? 7 r, 7 ii into the 
definitions of t,x,z,w, and thus we can set rj^ = fjr = 1 , without any loss of physical generality, which completes the 
proof. Furthermore, thanks to a suitable application of additional virtual symmetries S', we can choose the original 
phases r]r,riu such that = 1, where rjt = argt is the phase of t in (47): in this case we find that, without 

any loss of generality, we can set t > 0 , and even t > 0 (since for t = 0 the physical and virtual modes are decoupled). □ 


Next, we wish to incorporate translational invariance as well, on top of the other symmetries, leading to the final 
parametrization of the PEPS. Now we have translational invariance with unit cell size of 2 x 2, but we wish to 
incorporate an invariance without blocking, such that translational invariance also means charge conjugation: we 
require that the state will be invariant to first translating the particles by a single lattice site (which is the fermionic 
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part of the charge conjugation) and afterwards exchanging the positive and negative virtual modes, which is an 
inversion of the virtual electric field (to become physical once we introduce bosons and make the gauge symmetry 
local). 

For that, we decompose the operator A for even and odd vertices as follows: 


^(x) 


(x)gr„ a+ (x)6t (x) ^ 
(x)a] (x) gTij bt (x)at (x) ^ 


(48) 


We now act with the operator Ut (ei) on the state \ip (T)). The result is a somewhat peculiar state, in which the 
virtual fermions at x are connected with physical fermions at x + ei: 


Ut (ei) lip) = (12^1 (x)?^(x)A(xp = x + ei,x„ = x) |f2) 

X 


(49) 


We also define a similar translation operator for the virtual fermions, C4,. If 






(50) 


we obtain that acting with Ut (ei) on the physical level is equivalent to acting with C /1 (ei) on the virtual level, 
and similarly for the other direction 62 ; then, since the projectors and the vacuum are invariant under the virtual 
translation, we deduce that if (50) is fulfilled, the state \ip) is translationally invariant if T has the form 


T = 


/ 

t 


Ppt 

Vpt 


0 

y 

zl^/2 

z/V2 


-y 

0 

-z/y/2 

z/V2 


-z/V2 

z/V2 

0 

y 

V- 

-z/V2 

-z/V2 

-y 

0 


(51) 


where z,y &£, and t > 0, rjp = and the different normalization of 2 ; is chosen for later convenience. 

No further phases could be reduced, since we have exploited all the virtual redundancy symmetries which respect 
translational and rotational invariance. Thus, we arrive at the final parametrization of our PEPS, which shall be 
given in the following summarizing statement. 


Statement 3. The most general fermionic Gaussian PEPS with two virtual fermions per bond, with global U{1) 
gauge invariance, translational invariance and rotational invariance may be parameterized by three parameters - t > 
0,y,z G C, and is given by (51). 


6. The Covariance Matrix 


After having characterized the state in terms of the matrix T in (51), we start investigating some of the main 


properties of the PEPS |'0(r)). In particular we will define, in the following, its covariance matrix, parent Hamiltonian 
and correlation functions. To this purpose, we define a set of Majorana operators associated to the physical fermions 
as: 


Cx.l = V'x + ; Cx,2 = * (^/’x - V'i) > 

and their Fourier transformed (complex) operators: 


(52) 


rfk,a = (53) 

X 

where we consider a system oi L x L lattice vertices with periodic boundary conditions. Note that for k = 
(0, 0), (tt, 0), (0, tt) , (tt, tt), the resulting operators are still Majorana operators, whereas for the rest they are complex 
fermions. The redundancy appears as d^ ^ = d-^^a- Since the state \'ip) is Gaussian |53j, it is fully described by its 
covariance matrix, 

r^r' = ^([Cx,a,Cxq6]) 


(54) 
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which may be calculated using a Gaussian mapping (see|^. This is a real, antisymmetric matrix fulhlling T^T < 1, 
where the equality occurs for pure states. As the state is translationally invariant, the covariance matrix may be 
decomposed into blocks in momentum space m, which we denote as 


(Gout (k)U=^( 


rl 

“k,aJ “k,b 


iP (k) Q (k) 
-Q(k) -zP(k) 


(55) 


where P (k) is a real function whereas Q (k) is a complex function which may be decomposed into its real and 
imaginary parts Q (k) = R (k) + il (k). Due to the properties of pure fermionic Gaussian states [54] G^ut (k) = -1, 
therefore the functions P, R and I are normalized in such a way that P^(k) + /^(k) + P^(k) = 1. 

Inserting = d_k,a into Eq. ([55]) gives Gout(-k) = Gout(k), implying that 


P(-k) = -P(k), 
P(-k) =P(k), 
/(-k) = -/(k). 


(56) 


As shown inj^ the functions P(k), P(k) and /(k) are fractions of trigonometric polynomials, Po(k), Ro(k), /o(k) 
and P(k), which have maximum order 4 in ki and k2 individually, of the form 


p(k) = i;(k) = /(k) = 

^ ^ P(k) ’ ^ ’ V(k) ’ ^ ’ P(k) 


(57) 


where P(k) = P(— k) > 0. Hence, the relation s (|^ also hold for Po(k ), Ro{k) and /o(k) respectively. Due to the 
normalization of Gout (k) , they fulfill V (k) = \/RI (k) + Iq (k) + Pq (k). 


7. The Parent Hamiltonian 


Starting from the unnormalized covariance matrix 


0(k) = 


iPo(k) Po(k)+/o(k) 


-Rq (k) + Jo (k) -iPo (k) 

one can define a local Parent Hamiltonian [551156] for the state |'!/'(r)): 

H = -- ^ ' Qab (k) dk,adk_i, 


(58) 


(59) 


k,a,& 


with the dispersion relation V (k). By construction, the many-body state \ijj (T)) is a ground state of the Hamiltonian 
H. To obtain the Hamiltonian in real space, we define the Fourier transform of g (k). 


flab (x) = XI e*‘'’'flab (k) 


P2 


with the normalization chosen such that: 


H ^ ^ ^ Bab Cx,aCx^,b 

x,x' ,a,b 


(60) 


(61) 


or, in terms of the physical fermionic operators tjjx' 

H = X^o (x' - x) (V'iV'x' + V'I'V'x - <5x,x') + X (^0 (®2) 


x,x' 


where 


Ao (x) = Po (x) - iio (x). (63) 

Due to the fact that Po(k), Po(k) and /o(k) are trigonometric polynomials of maximum order 4 in ki and k2 
individually, their Fourier transforms can be non-vanishing only for \xi ^ — * 1 , 2 ! A 4, giving rise to a hnite hopping 
range of the above Parent Hamiltonian. 
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We observe that, in order to fulfill the global invariance, the Hamiltonian (62) has to commute with the staggered 
charge (14). This will be verified later by explicit calculation. Here we observe that the conservation of the staggered 
charge for a generic Hamiltonian of the kind (62) is verified if Ro{x' — x) = 0 for x and x' belonging to different 
sublattices and Ao(x' — x) = 0 for x and x' in the same sublattice. Furthermore, the global U(l) gauge symmetry 
implies the conservation of the fermionic parity in the system that is indeed evident in the Hamiltonian (62), which 
describes a system of spinless fermions subject to a p-wave pairing interaction, as we will discuss in Sec. 


8. The Correlation Functions 


The covariance matrix may also be used for the calculation of the correlation functions of the physical fermionic 
operators ^x- In particular, substituting'01 = (cx,i+jCx,2)/2 we can express them in terms of the real-space covariance 
matrix F defined in Eq. ( [54| ). We have 

4 (-*nT'+ nr' - nr' - ^nr') m 

or simply 

(V'iV'x') = ^ ((^x'x -Rn - x)) (65) 

and similarly 

= - ^A(x' -x) = (p(x) - if (x)) . (66) 


The first kind of correlation (65) 
as it would not preserve the 


IS 


global 


(where i?, A,P, / are defined similarly to g in (60), as Fourier transforms), 
expected to vanish if the two vertices are of different parities (sublattices), 
symmetry; Similarly, the correlation (66) must vanish in the other case. Let us show this explicitly. 

A gauge transformation of a Majorana operator is simply an orthogonal rotation matrix in SO{2), as it corresponds 
to a U(l) (phase) transformation of the fermionic operators. As the symmetry is global, there is a single parameter 
for these transformations everywhere; however, the rotation will have two different orientations, due to the staggering. 
And thus, in general, one may write the gauge transformation as 


pXX 

^ ab 


n 


O X pxx 
bb'^ a'b' ■ 


(no summation on x,x'). Denote 


f COS0 

I sin 0 


— sin 0 

COS0 


(67) 


( 68 ) 


as the transformation matrix of an even vertex, and thus corresponds to the transformation matrix of an odd 
vertex. Now consider an infinitesimal transformation (0 = e ^ 1) of a subblock of F involving two sites of the same 
sublattice (without loss of generality we take them to be even, otherwise one has to invert the sign). Then, 




o 


x'T 



-/(k) P(k) 
-P(k) /(k) 


e-*k(x-x') + o (e^) 


(69) 


j^xx p,g invariant, since the overall state cannot get a global phase under the transformation. As the second 

term must vanish for an even x — x' separation, we deduce that both P (k) and I (k) (as the imaginary part of Q (k)) 
must contain only odd harmonics. 

Similarly, if we consider the covariance matrix elements for two sites on different sublattices, we obtain 


J.X 


0X J^XX 0X _ J^X 


‘2c ^ 

k 




(70) 


and thus R (k) must contain only even harmonics. 

The last property we wish to discuss is the behavior of the functions P, P, / - either in real or momentum space - 
under rotations. 
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Statement 4. Under a rotation A, for \ijj (T)), 

R (/lx) = i? (x), R (ylk) = R (k) 

A (/lx) = —iA (x), A (/Ik) = —iA (k) 


Proof: From Eq. (65) we have: 


and, on the other hand, 


(V’axV’^x') = ^ (^yix',Ax -RiA (x' - x))) 


V’LV’Ax') = {Up {A)il)Uil)^,Ul (4)). 


(71) 


(72) 


(73) 


Since the state |^/> (T)) is invariant under rotations and satisfies Eq. (13) due to statement]^ it follows that 
(V'iV'x') which implies Eq. 

In an analogous way, concerning the correlator A we obtain form Eqs. (66) and ( fnpt : 


V’LV’ax') = -x)) = -*^t7p(/l)V'iV'I'C^p (7l)) ■ 


and Eq. ( [7T| ) follows again from the rotational invariance of |'0(F)). 

In particular, Eq. 0 implies that 

P(4k) = -/(k), /(4k) = P(k), 


(74) 


(75) 


and similar relations hold for the real-space functions. Furthermore, in accordance with the odd parity of these 
functions, we arrive at: 


P (-k) = P (/l^k) = -I (ylk) = -P (k) 
From the rotational invariance of |'0(T)) one obtains also 

d (/Ik) = d (k), 


(76) 

(77) 


as the dispersion relation should be rotationally invariant as well, and therefore: 


Po (^k) = Po (k), 

lo (vlk) = Po (k), (78) 

Po (4k) = -lo (k), 


along with similar relations in real space, which imply that the parent Hamiltonian (62) is also symmetric under 
rotations. □ 

The rotation relation linking P and / is a consequence of the simultaneous requirements of conservation of the 
staggered charge (14) and rotational invariance. In particular, as we will discuss in the next subsection, P and I 
correspond to hopping amplitudes of the fermionic matter after a staggered particle-hole transformation that maps 
the real space Hamiltonian (62) into a standard U(l) invariant Hamiltonian with a conserved number of fermions, 
reminiscent of the Kogut-Susskind model. The transformation rules of P and / under rotations are indeed consistent 
with the differences of the hopping amplitudes in the Kogut-Susskind Hamiltonian for staggered fermions along 
different directions (resulting from the different Dirac matrices) |44j . An example of these symmetry properties can 
be seen in figure 


B. The Globally Invariant PEPS as a p-wave Paired State: An Analytical Treatment 

1. The PEPS-BCS State in Momentum Space 

So far, we have considered some of the general properties of the correlation functions and the parent Hamiltonian 
by exploiting the Gaussian formalism and the physical Majorana modes. In this Section we calculate their explicit 
form as a function of the parameters in the T matrix. 
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2 

0 

-2 


R{k) 



-2 0 2 


2 

0 

-2 



-2 0 2 


2 


/(k) 



-2 0 2 


J(ylk) = P(k) 



-2 0 2 


P(k) 



-2 0 2 


P (ylk) = -I (k) 



-2 0 2 


FIG. 4. An example, with randomly generated parameters {t = 0.2785,1/ = 0.5469,2 = 0.9575), for the momentum-space 
matrix elements of the covariance matrix, showing explicitly the rotational symmetry. The colors represent the amplitude of 
the functions and are only important for illustrative purposes. N = 200. 


Before entering the PEPS details again, let us examine further some symmetry properties of the parent Hamiltonian 
in momentum space. For that, we adopt the following convention for the Fourier transform of the physical fermionic 
operators: 




\fLiL'. 


E' 


,ikx 




(79) 


where Li and L 2 are the 
particles. 

Using this convention, 
Hamiltonian: 


system width and length. Note that the momentum operators i/ik mix particles and anti- 


and the definition of Nambu spinors ifk = ('0k:V'Lk)^> obtain, from Eq. (62), the 


H = (k) (Po (k) a, + Jo (k) ay + Pq (k) a,) T-k. (80) 

k 


The Bogoliubov-de Gennes (BdG) Hamiltonian % (k) describes a spinless complex p-wave superconductor with order 
parameter Ao(k) = Pq (k) — Hq (k). A rotation {H —>■ Up (H) HU^ (A)) corresponds to 

n(k) ^w^n{k)w, with w = (si) 


This leads to 


W^n {k)W = n (Ak) 


which results in the rotation rules (78) discussed above. 

Such a p-wave pairing Hamiltonian enjoys the following particle-hole symmetry: 

ax'H{k)ax = -U (-k) 


(82) 


(83) 


whereas there is no time-reversal symmetry; therefore, in the general case, the Hamiltonian defines a system in the 
topological class D within the classification of topological insulators and superconductors (see, for example, [57H55] ). 

As J7 is a BdG Hamiltonian, its ground state, |i/> (T)) is a BGS state pairing modes with momenta k and — k, 
having the normalized form 

(u (k) -b V (k) i/ij,i/iik) l^^p) 

k>0 


(84) 
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2mk)|^ = (l-i?(k)), 

2u (k) z) (k) = A (k). 


These relations can be obtained by either diagonalizing the BdG Hamiltonian or with an analytical evaluation of the 
correlation functions, given by the covariance matrix obtained in the Gaussian mapping. In we follow the latter 
strategy to find explicit formulae for P, R and I as functions of the parameters of the matrix T. 

The PEPS \tjj{T)), however, is not normalized, and thus takes the form 


|z/>(k)) = (a(k)+/3(k)V'kV'-k) l^k) ( 86 ) 

where the functions a (k), /3 (k) being are unnormalized versions of the BGS functions u (k), v (k) and can be explicitly 
calculated with the Gaussian formalism. 

One may obtain from the previous equations: 


^ M (k) a (k) P (k) + il (k) ’ 


R{k) 

A(k) 


|a(k)|^-|^(k)|^ 

|a(k)|V|/3(k)r 

2d (k)/3(k) 


|a(k)|V|/3(k)r 


(87) 


where we defined the pairing function in momentum space g(k). 

Four modes, k = (0,0), (tt, tt) , (tt, 0), (0, tt), are left unpaired: these modes remain in the vacuum state (e.g. 
It/j (0, 0)) = d (0,0) |n)), i.e. jS = Q for these momenta. This means, in particular, that i? (k = 0 ) = 1 and A (k = 0 ) = 
0, in accordance with the previous results. The value of the pairing function g(k) in these points is instead non¬ 
trivial as we will discuss in the following. This particular behavior of the unpaired modes is dictated by the PEPS 
construction: occupying these unpaired modes is indeed impossible in |' 0 (T)), since the state was created by the 
operators A which involve only even products of creation operators. Despite that, we stress that, in the PEPS 
construction, a well-defined amplitude d is associated to these unpaired states (see[C|). 

Finally, the PEPS construction allows us to obtain a full analytical understanding of the globally invariant state, 
and to calculate explicitly all these functions. In particular the following results hold: 


Statement 5. The dispersion relation V (k) is proportional to 


E(k) = |a(k)|V|/3(k)^ 


( 88 ) 


which is the dispersion relation we shall work with. 
Statement 6. The BCS coefficients are 


a{k^, ky) = Aq + Ai [cos {k^ P ky) + cos {k^ - ky)] 

+ A 2 [cos { 2 kx) + cos { 2 ky)] + A 3 [cos { 2 kx + 2 ky) + cos { 2 kx 


2ky)] (89) 


and 


Pfkx^ ky) = 2t^Bi (sin kx — i sin ky) + 

2 t^B2 [sin {kx + 2ky) — i sin {ky — 2kx)] — 2 t^i ?2 [sin {2ky — fc^,) -f i sin {ky + 2kx)] (90) 

From these equations it is evident that: 


a{Ak.) = Q:(k), /3(2lk) = —z/3(k) 


(91) 
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The coefficients Ai and Bi depend only on y and z: 

Ao = (1 + + 3(1 + 2y^)z'^ - 4y^z^ + z® 

= -2 z2 (l + y^ + z^- 2y2 (1 + z^)) 

^2 = 4/z2 _ 2 y2 (^4 + 1) + ^4 _ 2^6 

A, = i - 2,^ f 

Bi = {l + z^){l + z^- 2yz^) + y^{l + 2z^ - z^) + 2zy^{2z^ - 1) - - 1) + y®(y - 2z) 

B 2 = y{z - y){l +y^- Z^) + ^z{-2y + 2y® + z - 3y^z + z®) 

B; = y{z - y){l + 2/2 _ ^2) - y{-2y + 2//® + z - iy'^z + z®) 


(92) 


Statement 7. The amplitude a associated to the unpaired modes |'(/'(ko)) = a(ko)|ri) of the PEPS at the momenta 
ko = ( 0 , 0 ), (tt, tt) , (tt, 0 ), ( 0 , tt) is related to the coefficient a of the paired modes by 


a (ko) = \/a (ko). (93) 

The proofs of these three statements rely on the PEPS construction, and are given in[C| 

We also observe, from the previous results, that a (k) >0 for all the values of the parameters y and z, and its 
minima lie at either k = (0,0), (tt, tt) or k = (0, tt), (tt, 0). In particular, the extreme values of a are 


oo = a(0,0) = a(7r, tt) = a‘^{0, 0) = (l - (// + zff {l - {y - zf^ > 
= a(7r, 0) = a(0, tt) = a‘^{Tr, 0) = (1 — //^ + z'^)^ 


Following these results, using eqs. (87), we can redefine the Hamiltonian coefficients (up to a constant) as 


i?o(k) = |a(k)|2 - |/3(k)|2 , Ao(k) = 2a(k)/3(k) 


in full agreement with the spectrum of the BdG Hamiltonian 

E (k) = ±^i?2(k) + p2(k) + /2(k) = ± (|a(k)|2 + |/3(k)|2) . 


(94) 


(95) 


2. The Phase Diagram 


Exploiting the analytical expressions (89) and ( |90[ ), we present here the detailed analysis of the phase diagram 
associated with the globally invariant PEPS \'if {T)), which is shown in Fig. Such an analysis involves, in principle, 
three parameters, f > 0 and y, z € C. t, however, is irrelevant in the definition of the thermodynamical phases: the 
closing of the gap E(k) may occur only at the four momenta defining the unpaired states, and it does not depend on 
t since a (k) is independent of this parameter and (3 (k), which is proportional to f^, is always zero in these points of 
the Brillouin zone. Therefore only y and z are relevant. 

For the sake of simplicity, we will mainly refer to the phase diagram in the plane defined by ?/, z; € R. We wish to 
emphasize, though, that the previous definitions of both the functions a (k) and /? (k) and the spectrum E (k) are 
valid for any //, z G C. Thus, the extension of the phase diagram to complex values of y and z is straightforward. 

As argued above, band touching points may only exist for the unpaired momenta, k = (0, 0), (tt, tt) , (tt, 0), (0, tt). 
The system becomes gapless {E (k) = 0 for some k) if and only if oq = 0 or = 0 [see Eqs. ([^]. Oq = 0 along 
the four lines z = ±1 ± y, whereas aT^ — 0 along the two branches of the hyperbola y^ — z^ = 1. In the first case 
band touching points will appear at k = (0, 0), (tt, tt), in the second for k = (0, tt), (tt, 0). These results hold also for 
complex values of y and z. 

To investigate further the characteristics of the system, we observe that the functions u (k) and v (k) are related to 
the Hamiltonian by 


Kk)|2 = ifi + ^' 




E{k) J ’ 




(96) 


Consistently with the usual analysis of p-wave superconducting systems, we determine the phases appearing in the 
model as a function of y and z, by considering the gap and the behavior of u (k) and v (k) in its minima. The minima 
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FIG. 5. The phase diagram of the fermionic, globally invariant model. The lines represent the gapless phases (A)-(E) described 
in the text, while the other areas are gapped. As explained in the text, this is only the plane oiy,z £ R, but the results are 
valid for an analytical continuation for any y, 2 G C. The Chern numbers of the gapless lines are indicated along them. They 
were calculated from the vector (Pq (k) , 7 o (k), Po (k)) characteruzubg the BdG Hamiltonian (8O1. 


of the dispersion relation are given by oq or therefore these parameters are related to the chemical potential in 
the system. In particular, since a (k) > 0, the chemical potential is always negative. 


Gapped Regions of the Phase Diagram 


In all gapped regions, where z ±y ± 1 and 1, R — Rq/E —>■ 1 for k in the minima of the energy 

iil(k). This implies that in the gapped phases, in a neighborhood of these minima, u —>■ 1 and u —>■ 0. Following Read 
and Green |60], this corresponds to a gapped regime with a strong p-wave pairing, such that the Cooper pairs are 
localized. The opposite situation, with u —> 1 and u —>■ 0 can never be realized in the gapped phases of our model 
since oq (k), a^r (k) > 0. This implies also that non-trivial gapped topological phases cannot appear in the system, 
consistently with the previous studies of chiral fermionic PEPS [ssmi]. Indeed, the winding number of the spinor 
(^(k), ?;(k))‘'' in the Brillouin zone is always 0 in the gapped phases, because u (k) 0 everywhere. 

In all gapped phases, ^(k) in Eq. (871 is an analytic function of ki,k 2 , since a(k) > 0. This implies that in real 
space the pairing function y(x) has to decay faster than any inverse polynomial of |x| (see Theorem 3.2.9. in Ref. |Ii2)l. 

The general study for arbitrary y and z is nevertheless complicated. Therefore, we focus here on two specific cases, 
as paradigmatic examples of the gapped phases. The first is given by y = I and z = v^, which we call the magic 
point: In this point (or in the equivalent ones obtained for y = ±1 and z = ±y/2), a (k) becomes independent of the 
momentum, and one obtains: 


a (k) = 16 , 

(3 (k) = 16 ^2 — t^ (sin ki — i sin ^ 2 ) . 

By rescaling a (k) and /3 (k) by the uninfluential factor 16, one obtains 

i?o (k) =1 — ^2 — V 2 J (cos 2 ki + cos 2 / 02 ) , 

Pq (k) =2 ^2 — sin ki , 

Iq (k) =2 ^2 — sin /c 2 . 

and for the pairing function 


(97) 


y(k) oc (sin/ci — IsinA: 2 ), 


(98) 
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FIG. 6. Correlation functions in coordinate space (logarithmic plots) for the magic point y = 1, z = \/2, with = 50 and a 
randomly generated t = 0.6324. Local gauge invariance is manifested by the vanishing of the relevant correlation functions for 
odd/even separations x; Locality is manifested by the exponential decay, which is apparent taking the color scale into account. 


such that g(xi — X 2 ) is nonzero only when xi and X 2 are nearest neighbors. Therefore, only maximally local Cooper 
pairs appear. In Fig. one may see plots of the correlation functions in coordinate-space, which manifest both 
locality and global invariance, for the case of the magic point. 

When we consider different values of y and 2 ; in the gapped phases, we expect the superconducting pairs to be no 
longer maximally localized, but their space distribution becomes exponentially decaying in the distance between the 
two fermions, consistently with a strong-coupling. 

One may also easily consider a second case, the “trivial” one, in which y = z = 0, and, in fact, half of the virtual 
modes do not participate (since the virtual subblock r of the T matrix vanishes, only the virtual modes whose sign 
is the opposite of the physical modes participate). In this case, a = 1 and (3 = 2t^(sinfci — i sin k 2 ), therefore 

i?o (k) = I — + 2t^ (cos 2ki + cos 2 ^ 2 ) , 

Iq (k) = sin k2 , (99) 

Pq (k) = sin ki 


and the spectrum is 


i? (k) = ± (1 -I- — 2t'^ (cos(2fci) -I- cos(2A:2))) (100) 

which, as expected, is always gapped. 

In real space the corresponding p — ip pairing Hamiltonian reads 

H={1- 4t^) + H-C-) + + H-C-) (101) 

XX X 


If one applies the particle-hole transformation ^1 —>■ ij}^ on the odd sites, the previous Hamiltonian is transformed 
into: 

H={1- 4f^) ^ (-1)“!+^= + H.C.) 

X X 

- 2f ^ (^iV'iV'x+ei + (-1)^'+"^" V'iV’x+es + H.C.) . (102) 

X 

The first and last terms can easily be recognized as the staggered mass and hopping terms of the 2-1-1 dimensional 
Hamiltonian for free Kogut-Susskind fermions [44j . In particular, we observe that the complex p-wave pairing gives 
rise to hopping phases that are consistent with the rotational requirements of the Kogut-Susskind model and constitute 
the remnants of the Dirac matrices in the continuum limit. 

The second term, instead, does not belong to this well-known model; furthermore we observe that these next- 
nearest-neighbor tunnelings cannot be minimally coupled to a static gauge field on the links of the lattice. 
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The spectrum coincides, non-surprisingly, with the one previously calculated, which means that the ground state is 
a state at half filling. For t = 0, all the particles occupy the odd sites (Dirac sea), corresponding to the vacuum state 
in the Kogut-Susskind Hamiltonian |44j . By increasing t, however, the even sites become progressively more and 
more populated. This corresponds to the formation of maximally localized pairs. Due to the form of T, indeed, the 
pairs cannot spread on the lattice, but they are composed of pairs of neighboring particles, in the superconducting 
picture, or by the displacement of single fermions from an odd to a neighboring even site for the number-conserving 
Hamiltonian (102). 


Gapless Regions of the Phase Diagram 

Let us start with a Taylor expansion of a (k) and /3 (k) around k —?> 0, 

/3 (k) =^i(fci - ik 2 ) + /33*ifc2(fci + ik 2 ) + Psikf - ikl) + 0 {k^) 
a (k) =ao + 02 {kj + fc|) -I- 0(A:^), 


(103) 


with 

^, = 2t^l + iy + zr) {1 -iy-zrf 

133 = -2it^[4y^ - 2y^z + z^ + z-^ + y^{4: - 7z^) + yi-6z + Az^)] 

/3' = I [2y(-9 + 8 y^ + y^)z + (7 - 28y^ + y^)z'^ + {y^ - 1)(1 - y^ - Ayz^) + (7 + y'^)z* + z^2y - z)] 


(104) 


a3 = {l-iy + zrr{l-{y-zrY (,05) 

a 2 = 2 [ 2 y^ + z\-l + z^f + 2 y\l + z^)-y\A + 3z^)]. 

For ao ^ 0 we recover that 5 (k) oc (fci — *^ 2 ), therefore the system is in a strongly paired gapped phase with 
the Cooper pair wavefunction, proportional to g(x), decaying exponentially in real space m- We distinguish the 
behavior of the gapless lines according to the values of the ai,fii coefficients, and identify five different regions: 

1. Strong Pairing Gapless Lines (A) : z = y 1, except for the points (y = 0, z = ±1) and {y = ±l,z = 0). 
Here both ag and /3i disappear. In this case the leading term of (/(k —>• 0) is of order 0(fc), corresponding, 
at large distances, to a real space pairing function dominated by |g(x)| = l/|xp. Since the average distance 
between paired fermions is finite, this is again consistent with a picture of localized Cooper pairs, thus with a 
strong-coupling regime. Furthermore, in this case too, u (k) —> 1 for k —> 0, and the Chern number associated 
to the spinor is trivial. 

The dispersion relation for k —>■ 0 is expanded, in this case, as 

£;(k)«a?fc4 (106) 

- a quartic dispersion in terms of fc = |k|. 

2. Weak Pairing Gapless Lines (B): z = —y ± 1, except for the intersection points {y = 0, z = ±1) and {y = 
±1, z = 0). Here, oq — 0,/3i ^ 0. In this case the leading behavior of g (k) is proportional to l/(fci -I- jfe), since 
the linear term /3i does not vanish. Therefore g(x) oc (xi +ix 2 )~^ for large distances and the system is, indeed, 
in a weak pairing p-wave gapless phase. In particular, since /3 dominates over a for small momenta, we obtain 
V (k) —>■ 1 and u (k) —)■ 0, and the Chern number becomes non-trivial (-2). 

In this phase, the behavior of the band touching points is described by 

E (k) « ^Ik^ (107) 

- a quadratic dispersion relation. Furthermore, close to k = 0, the leading behavior of the superconducting 
order parameter is given by 

Aq oc a2(3ik^{ki — ik 2 ) (108) 


and therefore the dominant pairing is of the pi — ip 2 type. 
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3. The Gapless branches of the hyperbola = 1 (C), except for the points {y = ±1, z = 0). For = 1 

the band touching points are k = (0, tt), ( 7 r, 0 ). Considering the expressions for /3 in Eq. ( [^ , it is easy to see 
that a translation of ^2 (or, analogously kf) by tt in momentum space corresponds to a transformation from a 
Px — ipy superconductor to a Px + ipy superconductor. Such a translation simply corresponds to the mapping 
ipx —>■ (—which does not affect |g(x)| in real space. 

Let us consider the band touching point at (0, tt) for — z^ = 1 and expand a and /3 in series of the momenta 
around that point (we define ky = ky — n): 


a (k) « z^{kf + k^) + 2(8 + 8 z^ + z'^)k\k 2 > 

ft (k) « —4it^{yz — z^ — l)(z^ + 4 )kik2{ki — ik 2 ) + 4t^z^(l — yz + z^) (jtf + ik^ , 


(109) 


where the first terms of the series vanish due to y^ — z^ = 1. With the exception of the points y = ±1, z = 0, the 
pairing function g()i) around the band-touching points has always a non-analytical behavior of the kind 1 /fc, 
therefore these critical hyperbola branches describe weak pairing gapless phases. 

The leading behavior of the dispersion relation close to the band touching points k = (0, tt), (tt, 0) is of order 
fc®. The Chern number along the hyperbola is nontrivial again (-1-2). 

4. Weak Coupling Intersection Points, {y = 0, z = ±1), (D). Here, Uq = a 2 = fti = 0, and therefore the leading 
term of g(k) is of order 0{l/k). This brings the system into a weak pairing gapless phase where the wavefunctions 
of the Cooper pairs decay as the inverse of the distance in real space. The Chern number here is calculated to 
be - 2 . 

5. Weak Coupling Intersection Points, {y = ±l,z = 0), (E). These are the intersections between the critical 

straight lines and hyperbola branches. These are the only points for which the gap closes at both k = (0, 0) 
and k = ( 0 , 7 r) since ao = = 0. We focus on the case {y = I,z = 0): The behavior of the other point, 

[y = —1, z = 0 ) is similar. 

The functions a (k) and ft (k) are, in this case. 


a (k) = 16sin^ ki sin^ k 2 , 

ft (k) = sin ki sin A: 2 (sin k 2 — i sin fci), 


Due to their common factor IGsin/ci sinA: 2 ) the Hamiltonian assumes the form 

% (k) = 64 (sin^ ki sin^ k 2 ) [Rq (k) Gx + /q (k) <Jy + Pq (k) = (sin^ fci sin^ ^ 2 ) % (k) 


( 110 ) 


( 111 ) 


with 


Po (k) = 1 — 4^"^ -I- (2t^ — 1) (cos 2ki + cos 2 k 2 ) -f - [cos (2fci -|- 2 k 2 ) + cos (2fci — 2 fc 2 )] 

/q (k) = 4 t^sinfc 2 — 2t^ [sin(2A:i -I- /c 2 ) -I-sin (^2 — 2ki)\ , 
pQ (k) = AP sin ki — 2 p [sin ( 2 k 2 -I- fci) -|- sin (fci — 2 ^ 2 )] 


( 112 ) 


The spectrum of TL (k) is 

if = ± [1 -I- 4t^ — (1 -I- 2P) (cos(2fci) -I- cos( 2 A: 2 )) + cos(2fci) cos( 2 fc 2 )] 


(113) 


which is always gapless with quadratic band touching points in both k = (0, 0) and k = (0, tt) for every value of 
t. Therefore, the leading term of the spectrum of TL (k) for k —>■ 0 satisfies 


E « 256Pkik2{ki + k^) 

and an analogous behavior is found close to (0, tt). 

Let us now evaluate the pairing function: 

(sinfc 2 — isinfci) P 


(114) 


5(k) = 


sin fci 


— i 


sin fco 


sin fci sin fc 2 


(115) 
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This corresponds to a weak pairing gapless phase (independent of the considered minima). We observe that, since 
5 (k) splits into 51 (fci) + 32 (^ 2 ) we obtain, as expected in the case z = 0, that g(x.) = gi{xi)Sx 2 ,o + 52 (a^ 2 )^xi,o; in 
particular g(x) oc —it'^{—l)^^Sx 2 ,o — for x ^ (0,0). Therefore the Cooper pairs, which correspond 

to the mesons in the gauge theory, do not decay with the distance between particle and antiparticle but are 
constrained to spread only in the horizontal or vertical direction. The Chern number at these points is trivial. 

The phase diagram is summarized in Fig. 


IV. LOCAL GAUGE SYMMETRY 


So far, we have introduced purely fermionic states \tp (T)) which, on top of the spacetime symmetries corresponding 
to translation and rotation invariance, are invariant under global U{1) transformations. Now we wish to lift the global 
symmetry to be local, i.e. instead of the global transformation rule with a global transformation parameter 

(phase) (j), we would like to have a local transformation rule, with vertex-dependent parameters (phases) ())x, 

^ (116) 


where Sx = (—l)“i +®2 accounts for the staggered charge of the fermions. In general, the states \ip {T)) defined in 
the previous section will not be invariant under such a local gauge transformation and need to be modified. This 
modihcation is necessary because the effect of the transformation of the phase of a single vertex ( |116 ) must be 
compensated by the action on some additional degrees of freedom which cannot depend on other vertices if the 
symmetry is local. Geometry and locality considerations lead us to introduce, as customary in lattice gauge theories, 
new degrees of freedom on the links connecting the vertices, which participate in the local gauge transformations in 
a way that conserves the local charge (which will be the eigenstate of a physical Gauss law). 

Motivated by compact QED (cQED) [7], we introduce a bosonic Hilbert space on each link, and follow a procedure 
similar to [39]. These Hilbert spaces shall be denoted as "H® and "Hx, representing the links on the right (side) and 
above (top) the vertex x (see Figure]^. As we wish these Hilbert spaces to be finite, we shall work with those appearing 
in truncated cQED linillTj: the Hilbert spaces are spanned by the SU{2) states {|w)}^^_^ = for 

a given, fixed integer £ (full cQED is obtained for £ —>• 00 |T7]1. 

These states are eigenstates of the electric field, T, = L^. We define the raising and lowering operators 


S± 


L± 

vW+i) 


(117) 


having the t/(l) limit 


|m) —|m) = |m -I- 1) (118) 

^—>•00 

where 0 is the cQED “vector potential” conjugate to the electric field. 

Therefore, for each lattice vertex x, we consider the Hilbert space of the states of both the matter fermion sitting 
on the vertex and the two gauge bosons on its right and top links. Within this extended space, we have to consider 
the symmetry requirements associated to the new, locally gauge invariant states, \'fb)i which will generalize the 
symmetries fulfilled by the purely fermionic states \ip). 


1. Charge Conjugation Symmetry. This is a generalization of translational invariance due to the effect of the 
staggered charge of the fermions. A translation from one sublattice to the other corresponds to a charge 
conjugation that has to affect the gauge bosons as well. Therefore we extend the action of the Ut operators 
defined for the fermions ([^ to include the bosonic operators acting on 77® and 77* as well 

Ut (Sfe) Sf (x) 4 (efc) = (x + e^), 

Ut (bfe) (x) 4 (efc) = -S®/* (x + e^). 

Such a definition of Ut introduces, in fact, a charge conjugation: it both exchanges matter particles and anti¬ 
particles, and inverts the sign of the electric field. We wish it to be a symmetry, and thus we demand 


Ut (efc) \fb {{U})) = |'0b ({£*})) • 


(120) 







FIG. 7. The Hilbert space on a vertex (composing the fiducial state), in the case of a locally gauge invariant state: a single 
physical fermion tp and two physical bosonic states, s and t, with eight virtual fermions surrounding it, two on each edge 
intersecting at the vertex. 


2. Rotational Invariance. The lattice rotation has to affect both the fermionic vertices and the bosonic links. We 


extend the fermionic rotation ( 11 ) with: 


where 


and the symmetry requirement reads 


C/p (H) Si (x)C/t(T) = Ei (Tx) 
C/p (H) Si (x) C/pt (T) = Si (Tx) 
C/p {A) E"/‘ (x) Ul (A) = E®/* (Tx) 


Si(x) = Si (x-ei) 
Si(x) = Ei (x-e2) 


Up (/I) IV'fa ({/*})) = IV'h ({/*})) • 


( 121 ) 


( 122 ) 


(123) 


3. Local Gauge Invariance. We finally define the local C/(l) gauge transformations, generated by the Gauss law 
operators, 


Gx = E^ (x) + E* (x) - E® (x) - E‘ (x) - 

= E^ (x) + E* (x) - E^ (x - ei) - E* (x - 62 ) - Q^. 


(124) 


This generator extends Eq. (18 1 to the gauge degrees of freedom. As one can see, the virtual electric fields Ei 


from the previous section have been lifted to physical degrees of freedom, and the virtual local Gauss’s laws 
which defined the fermionic fiducial states with a global physical symmetry, have turned into physical local 
conservation laws. Therefore the physical symmetry we demand for the locally gauge invariant state is 




I'^b {{U})) = \lpb {{ti})) Vx 


(125) 


Hereafter we focus on the case C = 1 which constitutes the simplest truncation scheme allowing for zero, positive 
and negative electric fluxes on the links. In this case the electric field states are simply |0) , |±1), and 

/O 1 0\ /O 0 0\ 

E+ = 0 0 1 ; E_ = 1 0 0 (126) 

\ 0 0 0 J \0 1 0 J 

in the basis {|+ 1 ), | 0 ), | — 1 )}. 

In the following, we will present a recipe to transform the globally-gauge-invariant PEPS of section in into a 
locally-gauge-invariant one, and show that the symmetries are obtained using the same parametrization presented 
for the purely fermionic case: one can parametrize \tjjb) with the same matrix T {t,y, z) introduced for the fermionic 
Gaussian state \4>). However, the bosonic state will not be Gaussian any longer, reflecting the physics of the underlying 
interacting theory. The reader could also refer to other PEPS constructions, as in for pure-gauge theories, or [3S] 
for gauge theories coupled to bosonic (Higgs) fields. 
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A. Inclusion of Gauge Bosons in the PEPS 

Following the steps of the construction of the fermionic I'i/'), we start by imposing the local gauge invariance for the 
fiducial state. On each vertex we define a fiducial state 

(x)) = Ab (x) |Op (x)) |0s (x)) |0t (x)) |0„ (x)) = Ab (x) |0 (x)). (127) 

with the operator Ab yet to be dehned. This will be a state involving one physical fermion at the vertex, and two 
bosonic degrees of freedom corresponding to the right and top links connected with the vertex, as shown in Fig. 
The corresponding PEPS takes the form: 


IV'b) = (x) 77 (x) Ab (x) | 0 ) 


(128) 


This physical state is analogous to the fermionic PEPS \ip) in (FtI), with the same projectors w (x), 77 (x), defined by 

©• 


Statement 8. The state \ijjb) is locally gauge invariant, i.e. it fulfills [125), if the following conditions on the fiducial 
state are met: 

1. U^U\r \Fb) = UtUsUl \Fb) 

II. U^y \Fb) = Us \Fb) 

III. \Fb) = Ut \Fb) 

where U., 1 , = and the virtual operators are defined the same as for the fermionic theory, whereas 

Us = and Ut = act on the physical bosons. The previous relations can be rephrased for the Ab operators as: 

t. U^UlrAbU\}U^^ = UtUsUlAbU^WsU} 

ii. UifAbUD = UsAbWs 

III. U^AbU^^ = UtAbUl 

Proof: Local gauge invariance, written with the operators of the statement, is 

rt rrt T-rt u;,. \ _ tt tt. tA tA tA 


IV'b) = C/.,xC/7,xC/I,x-e,C^/.x-e.t^;.x l^b) = IV®b) 


(129) 


for any x. Using the PEPS definition (128), we find out that if the properties of the statement are fulfilled, then 


IV'b) = (y) V (y) Us (x) Ut (x) u] (x - ei) ul (x - 62) ul (x) \Fb (y)) 

y 

= (U„| (y) ?7 (y) U^ (x) U^ (x) Uy^ (x - ei) (x - @ 2 ) (x) IF^ (y)) = li/jb) 


(130) 


where the second equality is due to to statement’s conditions, and the third one due to the local invariance of the 
projectors under the transformation (28). □ 

In order to guarantee the local gauge invariance of \'ipb), we have to construct a fiducial state fulfilling the conditions 
(I-III). To this purpose, we will adopt the following construction: we consider the fermionic state \ip (T)), and lift the 
virtual degrees of freedom to be physical in order to enforce the local symmetry defined by the physical generator Gx 
in (124), by, as we have mentioned, lifting the virtual electric fluxes of Go (as in ([l4|)) into physical ones. Eor that, 
to define Ab, we multiply the virtual fermionic operators of the right and up edges by physical bosonic operators as 
follows: 


rl 

u[_ —)■ 'E^_u[_ 


(131) 
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Let us verify that the properties of statement are fulfilled. Properties ii and iii are easily verified: Since both 
UyT^JJy = and one obtains that Uyrl.T,^Uy = Usr^S^C/J' and Property ii is deduced. 

Property iii may be similarly shown. 

The first property requires more care. Using the Gaussian construction of the fermionic part, accord ing to statement 
we have \Fh) = \Fi,), where Go is the local fermionic virtual Gauss operator defined in (18). Therefore 
UyUy |Fb) = UyUyll^ jU^). Usmg propertfes ii and iii of statement which have already been proven, one obtains 
property i. 


1. Rotational Invariance and Charge Conjugation 


The next symmetry we check is the rotational invariance. 

Statement 9. The states |'0t, (T)), obtained with the same parametrization as of {ip (T)) by making the replacements 
(131), are rotationally invariant. 


Proof: Combining the fermionic transformation rules (11),(36) with the bosonic ones (121), one gets 

Up (T) Ur (vl) Ab (x) Ul (T) ul (^) = A( (Tx) 


(132) 


where is an operator involving a bosonic Hilbert space attached to the left edge rather than to the right, i.e. with 
physical states t and s. 

However, thanks to the horizontal projectors w, we can convert into Ai,: 


{a 


uRow I 


a;(x) (x) |U„Row) — (U^Rowl UJ (x) n A.I) (x) |rii;Row) 


xGRow xGRow xGRow xGRow 

and similarly in the vertical direction. Thus we may conclude that, once T is parameterized properly, 

Up{A)\(;, (r)) = |v>6 (T)) 


(133) 


(134) 


which proves the statement. □ 

The last remaining symmetry is charge conjugation. We have already defined this transformation as a generalization 
of the fermionic translation. Combining the fermionic translation (|^ with the bosonic electric field inversion ( |119 ) 
may be straightforwardly understood as charge conjugation, i.e., both exchanging particles and anti-particles (by the 
translation) and the sign of the electric field. The fermionic part is guaranteed using the parametrization introduced 
before. All these results lead us to the final statement. 


Statement 10. The state \^pbiT)), defined by making the replacements (131) in the fermionic state {ip (T)) of state¬ 
ment^ having the same parametrization matrix T = T (t,y,z) with t > 0,y, z G C as in Eq. (51), has a local U{1) 
gauge invariance, as well as rotational invariance and charge conjugation symmetry. 


B. The transfer matrix associated with the PEPS 


For the numerical implementation of the locally-gauge invariant PEPS we adopt a cylindrical geometry, with 
periodic boundary conditions in the ei direction and open boundary conditions in the 62 direction, consistent with 
having all the bosonic states on the lower boundary set to jO)?,. We label by Li and L 2 the numbers of horizontal and 
vertical matter sites of the lattice, respectively. We follow the construction presented in for the study of chiral 
fermionic systems in order to study the PEPS \ipb) in (128) . In the following we present some of its main features. In 
particular, in our case, the presence of two virtual fermionic modes per bond implies a bond dimension 4. As already 
discussed, the state \tpb) is obtained starting from the product of all the fiducial states \Fb) on which the projectors 
u! and r] are applied. One may define, for the whole lattice, the state: 


n 


012 = 1 


n ) n 


|U(x)) 


\xi—l 


Xi —1 


Xi — 1 


(135) 


This state involves all the virtual and physical states in the lattice. 
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The expectation value of a physical observable Oy supported on the y-th row can be expressed as: 

{Oy)=tT[XjOyX,\T){T\] (136) 

where Xi = 1)| is simply a projector, acting on the first row of the virtual states associated to 

the modes d±{x 2 = 1), that fixes the boundary conditions in such a way that there is no ingoing virtual electric flux 
Ed in the Hrst row of the cylinder, which can also be interpreted as the additional presence of an initial row of gauge 
bosons in the state |0b), and is a similar projector setting the boundary condition for the upper row. 

To evaluate the previous expression, and in particular the density matrix it is convenient to introduce the 

reduced density matrix X{x 2 ) associated with the virtual d± modes of the row X 2 + ^- In this way, the evaluation 
of the expectation value Oy can be performed by dividing the calculation into one step per row. The density matrix 
X{x 2 ) can be defined iteratively starting from the corresponding density matrix in the previous row |63j . 


X{x 2 ) = tr 


W_'n{xi,X2)u:{xi,X2)Ab{xi,X2) |fl)(fl|(xi,a;2) 


Y\_Al{xi,X2)uj{xi,X2)ri{xi,X2) X{x2 - 1 ) 


(137) 


(note that the definition of trace in a fermionic Hilbert space is non-trivial, due to the nature of fermionic Hilbert 
spaces mentioned before, and thus has to be done with caution) where the projectors |H)(H|(xi, 0 : 2 ) are on the vacuum 
of all its physical and virtual fermionic modes and on the |0) bosonic states associated to the site {xi,X 2 )- The trace 
is taken over all the physical fermionic and bosonic states in the row X 2 , and over all the virtual fermionic modes, 
in such a way that X(x 2 ) is indeed an operator acting on the virtual modes d±{xi,X 2 + 1) through the projectors 
r]{xi,X 2 )- We have exploited the fact that, as projectors, w = and r] = rj^ (see [B5] for more details). 

The previous equation describes a mapping between X{x 2 — 1) and X{x 2 ) that can be expressed in terms of a 
transfer matrix T defining the mapping between the virtual density matrices: 


In particular we have 


x{x2)dA' = r['!x{x2 - i)u.. 


r = tr 


p,t,s,r,l,u 





(138) 


(139) 


where the trace is taken over all the modes associated to the row X 2 with the exception of the virtual d± states, in 
such a way that T acts on both the d virtual states of the rows X 2 and 0:2 + 1 . 

We can now rephrase the expression for the expectation value (136) of a local observable Oy in the following way: 


tr 


{Oy) = 


XfT^--yfojy-^x, 


tr [XfT^-X,] 


(140) 


where we have exploited that the observable Oy has support on the row y only. Toy is the specific transfer matrix 
associated to the row on which the operator Oy is acting. It is defined by introducing the physical observable Oy 
within the definition of T in Eq. (139) 


foy = tr 


Oy 





(141) 


This way of evaluating the expectation values emphasizes the role of the transfer matrix T which is fully defined by 
the operators Ah and can be numerically evaluated. The non-degeneracy of the maximal eigenvalue of T is related to 
the presence of a hnite correlation length in the physical system. 

It is believed that for PEPS a gap between the two highest eigenvalues of the transfer operator T indicates expo¬ 
nentially decaying correlations in real space. This has been shown rigorously for MPS [64] and is based on the fact 
that the correlation between any two local operators scales as the second highest eigenvalue to the power of their 
distance (the highest eigenvalue has to be normalized to I). This observation was also made in numerical studies with 
PEPS, though a rigorous proof is still awaited. 
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FIG. 8. Left: Gap between the highest eigenvalue Ai = 1 and the second highest eigenvalue of the transfer matrix for t = 0 
(top) and t = 1 (bottom) for a cylinder of circumference Li = 6 as a function of y,z € R. Right: Same plot for Li = 8. The 
lines with low values of the gap indeed seem to be gapless lines in the thermodynamic limit, as the corresponding values of the 
gap are significantly lower for Li = 8 than for Li = 6. Points were the Lanczos algorithm for calculating the eigenvalues of the 
transfer operator did not converge are marked in white. 


C. The phase diagram of the locally gauge-invariant state 


The opera tors Ab defined by the equations ( 21|51 1 with the introduction of the bosonic operators through the 
substitution (131) completely define the locally gauge invariant state \ipb) as a function of the three parameters t>0 
and y, z G C. In particular, the parameter t couples the physical matter fermions with the virtual and bosonic modes. 
For t = 0 no physical fermionic creation operators contribute to At, therefore I'l/'b) becomes a bosonic state of the 
kind |^/’{))iinks|^^p)vertices! non-trivial on the bosonic links in the lattice, whereas all the fermionic matter sites remain 
in their vacuum state \flp). This limit thus corresponds to the pure U{\) truncated lattice gauge theory defined in a 
locally gauge invariant sector without any static charges on the lattice vertices. 

Due to this feature, in the case t = 0 only bosonic operators and observables are meaningful and the theory can 
be restricted only to the bosonic Hilbert space. For all the values t > 0, instead, the matter fermions appear in the 
definition of |'0b)i and this sudden addition of fermions has to be related, at least intuitively, to some discontinuity 
or non-analyticity in the limit t —>■ 0. When we deal with the combination of dynamical fermions and gauge fields we 
shall address this issue. 

The calculation of the spectrum of the transfer matrix T presented above allows us to detect the presence of critical 
lines and surfaces in the phase diagram of the states z)). In Fig. we plot two sections of the phase diagram 

for t = 0 and t = 1 and real values of y and z. The figure displays the gap A between the two largest eigenvalues of 
the transfer matrix. In the following, the term ’’gap” will refer to this. 

In the pure gauge theory, t = 0, the signs of y and z play no role, and can be eliminated using phase transformations 
of the form (30), with phases ±7r. In the following we shall restrict our discussion to j/, z G R, and thus it will be 
sufficient to consider y, z > 0 only. In this case, four gapped phases can be easily spotted due to the appearance 
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of several critical lines, defined by Z\ = 0 compatibly with the resolution of our numerical calculations. To verify 
the reliability of the existence of such critical lines in the thermodynamical limit, we have repeated the numerical 
calculation of the gap A for different system sizes, on cylinders of increasing width, which seem to show that the gap 
A is closing in the thermodynamic limit along certain lines in the phase diagram. 

The analysis and characterization of the gapped phases has to be based on suitable order parameters. In the 
following, we will consider separately the pure lattice gauge theory (t = 0) and the full model with fermionic matter 
and gauge fields at t > 0. 


1. The pure gauge theory 


In the pure bosonic model at t = 0, the gauge invariant physical observables which can be exploited to characterize 
the states are of two different kinds: bosonic operators based on the local electric field S, which commutes with the 
Gauss law, and closed string operators built based on the Wilson lines I]±. The second class commutes with the 
local Gauss law only if suitable products of these operators enter in such strings based on the convention about the 
orientation of the links on the lattice, which is dehned by the signs in Gx in Eq. ( |124[ ). 

The local expectation value of the electric field vanishes due to the imposed charge conjugation symmetry. The 
operators with q S IR a generic parameter related to the elementary charge of the model, allow to define ’t Hooft 
loops on closed paths on the dual square lattice (see Fig. [^, 

Gx,(g) = (142) 

bea-D xGi? 


where I? is a closed region on the square lattice delimited by the physical links on the loop of the dual lattice dV (see 
Fig. 1^. The first product represents the application of the gauge transformation to all the bosonic sites b 

surrounding the closed region V. The signs Sb = ±1 are defined by considering the orientation of the loop on the 

dual lattice, as depicted in Fig. _ _ 

The equality between the two products in Eq. (142) is guaranteed by the Gauss law (124): the’t Hooft loop is related 
to both the total matter charge X)xgx> included in the region V and the electric field flux SbS(b) 

associated to its contour dV. In particular, in a pure gauge theory, the expectation value of the’t Hooft loop Gx>{q) 
is just an exponential of the static charges enclosed in the region V. In our PEPS construction for t = 0, though, no 
static charge appears; thus {G^iq)) = 1 for all the closed regions delimited by a contractible loop dV on the dual 
lattice. The situation is different when considering a non-contractible loop on the cylinder (see Fig. |^. In this case 
the’t Hooft loop becomes 


Li 

Gg^(g) = . (143) 

By applying the Gauss law, and considering the absence of static charges, it is easy to show that (G^^{q)) is 
independent of the coordinate X 2 and it simply represents the electrical flux flowing along the surface of the cylinder 

= ;^E‘(cci,cr2). (144) 

xi 


Such flux, independent of X 2 , is defined by the boundary conditions adopted in the first row and it is zero in our 
numerical calculations (unless stated otherwise). More in general, G^^{q) is independent of local perturbations: All 
the loops on the dual lattice surrounding the cylinder once are equivalent and measure the electric flux defined by the 
projector operator Xi in Eq. (1361. 

The second class of gauge-invariant bosonic string operators is the one of the Wilson loops which can be thought 
of as the path ordering of the exponentiated vector potential integral exp A^dx^) along a loop C on the lattice. 


Following Eq. (118), in our truncated lattice gauge model we apply the substitution exp J) 
Ey*(x). Therefore the (clockwise) Wilson loop reads 


x+e„ 


A^dx ^) = e 




Wc = 

bee 


(145) 


where the choice of the operators Ej- depends on the orientation of the loop C: For links oriented upward or rightward 
one has to apply E_|_, whereas for downward or leftward links, E_ (see Fig. [^. It is important to notice that, in our 
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FIG. 9. Wilson and ’t Hooft loops are illustrated in a cylindric system of width Li = 10. Circles and squares correspond to 
physical sites associated to the gauge bosons and the matter fermions respectively. The grey circles depict periodic boundary 
conditions in the horizontal direction. Wilson loops are illustrated as red lines on the lattice: both a non-contractible loop and 
a 2 X 3 closed loop are shown, ’t Hooft loops on the dual lattice are represented by dashed blue lines: both a non-contractible 
and a 3 X 2 closed loop encircling the area 7? are shown. For each loop we specified the operators acting on the sites along each 
edge consistently with counterclockwise loops. 


model, the Wilson loops W are not unitary operators. Furthermore they do not commute in general with each other. 
This is due to the relation [S+, S_] = E valid for our truncation with £ = 1. 

We observe that a vertical Wilson line crossing the whole cylinder from the first to the last bosonic row is a further 
gauge invariant operator of the system and it does not commute with the non-contractible ’t Hooft loops. Despite 
that, though, we have verihed that our state does not show any topological degeneracy. Indeed the transfer matrix T 
is block-diagonal with the blocks labeled by on the ket layer and on the bra layer. Hence, if a density matrix X 
has fixed fluxes on ket and bra, it preserves them individually. The largest eigenvalues of different sectors of T labelled 
hy ^E on ket and bra are, however, non-degenerate, differently from what happens in topologically ordered models 
such as the toric code IHH| . The dominating eigenvalue is the one associated with the = 0 sector. A rough 
finite size scaling indicates that the second largest at <Pe = ±1 for ket and bra has a magnitude of 7.9% as compared 
to the former. The absence of topological order in spite of the presence of two non-commuting gauge-invariant string 
operators can be ascribed to the non-unitarity of the Wilson loop operator W. 

We have verified that the correlation between non-contractible Wilson loops along horizontal lines at X 2 and 
decays exponentially in all the gapped phases appearing in the phase diagram at t = 0, as expected by the non¬ 
degeneracy of the transfer matrix 




(146) 


for a ;2 » X 2 , where we have dehned the non-contractible Wilson loops as W^^{x 2 ) = Jla;! ^+( 2 ^ 112 : 2 ) (see Fig. 10). 

Away from the critical regions, one of the main distinctive characters of the phases in a pure lattice gauge theory 
is the exponential decay of the expectation value of the Wilson loops as a function of its dimension in the limit of 
large loops. Such exponential decay is typically dictated by either the area A of the loop as or its perimeter 

V as . The former behavior signals a confinement of the static charges, whereas the latter would be compatible 

with a deconfined phase of the static charges |31[71E7] (although in a pure compact QED, without a truncation, the 
latter phase cannot exist [3 EH IM!). To evaluate this behavior, however, the calculations need to be performed in 
the thermodynamic limit of the system and for large enough loops in order to properly distinguish the two behaviors 
and avoid finite size effects. 

Our cylindrical systems are limited in the periodic dimension to a finite size of Li < 9. For the pure gauge case 
of t = 0 the staggering plays no role and we can indeed exploit also odd values of the width Li. When we consider 
























31 



FIG. 10. Semilogarithmic plot of the exponential decay of the correlation function of two non-contractible 
a function of their distance h [see Eq. (146l]. The data plotted correspond to I + Z 2 )) — 


Wilson loops as 
(W'^°(30))^| in 


a system of size 6 (Z 2 + 40). The correlation of the two Wilson loops decays exponentially in all the phases. Due to the fast 
exponential decay, for larger values of the distance I 2 the numerical errors becomes too large to obtain reliable data. 



FIG. 11. A schematic plot of the phase diagram for the pure gauge theory (t = 0), with y,z > 0 (straightforwardly generalizable 
to any y,z £ R as explained in the text). The A,B phases seem to confine static charges, while the C,D phases seem to be 
deconfined. 


rectangular Wilson loops W{li, I 2 ) with width li and length Z 2 , the maximal distance between the two vertical edges 
is at most Li/2 due to the periodic boundary conditions, and, in our case, we were limited by h < 5. This extension 
of the loops is sufficient to evaluate their asymptotic behavior only in the presence of a large enough gap A of the 
transfer matrix and this implies that our calculations are reliable only far enough from the critical lines in the upper 
phase diagrams in Fig|^ Therefore, in the following calculations, we mainly consider points located deep inside the 
bulk of the gapped phases. Furthermore we label the phases A, B, C, D as schematically shown in Fig. [g 

Keeping this limitation in mind, let us analyze the numerical results for Li = 8, 9 (the two sets of data are totally 
consistent with each other). We have considered loops of length I 2 < 20 embedded in a cylinder with L 2 = I 2 + 40 in 
such a way that the loops have a distance 20 from both the top and the bottom edges of the system. 

The phase D clearly presents a perimeter law decay of the Wilson loop for the values of the y, z that we probed 
(see Fig. [12)3). This is consistent with the PEPS structure obtained in the extreme limit y —>■ 00 , z = 0: an expansion 
of the fiducial state of a single lattice vertex around z = 0, 1/y = 0 results in a superposition of the bosonic vacuum 
(in zeroth order), and O {y~^) terms involving either horizontal or vertical flux loops crossing the vertex (since y is 
the T parameter responsible for horizontal and vertical coupling between links). Since no flux ’’corners” (z terms) 
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A = kk P = 2(/i+/2) P = 2(/i+/,) 

FIG. 12. The area/perimeter law behavior of the Wilson loop, in the phases B,C,D for the pure gauge case {t = 0). It is 
clearly seen that phase B is governed by an area law, while phases C, D respect a strict perimeter law. Data was computed for 
Li = 8. The two lines emanating from the main line of the area law in phase B belong to the data sets of Zi = 6, 7 and I 2 > 12,9 
respectively. These could be perimeter corrections due to finite size effects, as the widths of these loops are comparable with 
those of the system, Li. 


are there, for the first nontrivial order, the most significant contributions to the physical state would be arbitrarily 
long (as the system size) flux lines, as expected within a deconfining phase, showing a perimeter law behavior for the 
Wilson loops. 

The appearance of the perimeter decay of the Wilson loop can be explained from the properties of the PEPS in 
the limit z —)■ 0. In fact, the operator at t = 0 and z —)■ 0 becomes the product of two separate operators acting 
on the horizontal and vertical links respectively. We have indeed 


Ab{t = 0, z = 0,?/) = exp 




zLei) Oexp y (uLdLEi 


= Ah A,, 


(147) 


This implies that the PEPS is decomposed a product state of one-dimensional systems 


\^^{t = 0,z = 0)) = 

r]{xi,X2) Ay{xx,X2)\^{xi)) (g)(f^.(^2)in u{xi,X2)W_Ah{xi,X2)\^{x2)) 

X\ X2 X2 X2 X\ X\ 

= |ver)®'^i|hor)®-^= (148) 

where we splitted all the vacuum states in the product of states for the degrees of freedom aligned along rows and 
colums. |V'b(^ = 0,z = 0)) is therefore a product state of Li identical vertical ID states, labelled by |ver), and L 2 
identical horizontal ID states, labelled by |hor). When we consider a rectangular Wilson loop Wc, only four of these 
states are involved, one for each edge. In particular we decompose the Wilson loop into Wilson lines acting on the 
four ID systems: Wq = WbotlDieftlTtoplDright where these ID operators have a structure such that ITtop = M^bot 
Wieft = We obtain 


{^CihM)) = |(hor|Wbot|lior)|^ |(ver|tTieft|ver)|^ . (149) 

Such an expectation value vanishes exactly at z = 0, because the finite Wilson lines violate the local Gauss law in each 
ID system. In particular, in the decoupled z = 0 limit, due to the periodic boundary conditions, each ID horizontal 
state is a cat state of the form 

|hor) cx (1 + |0,0,0,...) + |1,1,1,...)+ |-1, -1, -1,...). (150) 

In the thermodynamic limit Li —> 00 , the first term prevails for every value oi y ^ ±1. At y = ±1, instead, the three 
terms have a comparable amplitude such that the critical line originating at y = 1 and z = 0 is characterized by a 
symmetry breaking with three possible degenerate states of the corresponding parent Hamiltonian, whereas for the 
gapped phases only the state with no electric field survives. In the vertical direction the states |ver) have a similar 
behavior but are constrained by the boundary conditions. 

For small values of z, all these decoupled ID states are perturbed by the introduction of domain walls with an 
amplitude proportional to z. Such domain walls describe indeed corners of the electric flux in the overall 2D state, 
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which becomes a set of weakly coupled ID systems. The Wilson loop Wc is an operator which can be decomposed into 
four of these corners, each one with amplitude z, in such a way that the 2D Gauss law is not violated. Therefore its 
expectation value decays as Cz‘^e~^ ( 2 I 1 + 2 I 2 ) |-]^g lowest order in z, where the appearance of the perimeter law is due 

to the fact that each ID state has a gapped transfer matrix away from the critical points, thus each ID expectation 
value decays with the length of the related Wilson line. Hence, for z ^ 1, the decomposition of the PEPS into 
one-dimensional states implies that no area contribution to the decay can be present at low orders in z. This explains 
the appearance of a perimeter law for the Wilson loop in both the phases D and C. 

Our numerical results confirm that the behavior of the gapped phase C is also characterized by a perimeter decay, 
despite being clearly separated by the phase D. In this central region of the phase diagram, though, the finite size 
effects due to the limited width of the loops seem to be more relevant. This implies that only the data at Li = 8, 9 
with li having a maximal value 4 are reliable, whereas data taken at Li = 6 present large deviation with respect 
to the perimeter law. Furthermore, also for Li = 8,9 we see minor deviations from the perimeter decay due to the 
different widths of the loops. The perimeter law behavior may be seen in Fig. [I^. 

Such behavior of the phases C and D is fully compatible with a deconfinement of the static charges: in fact, no 
signature of an area law decay appears in these phases. 

Let us analyze the Wilson loops in all the phases more quantitatively: To limit the effect of the perimeter contri¬ 
bution and evaluate the eventual decay as a function of the area (dictated by ka) in all the four phases, we evaluate 
the following parameter, introduced by Creutz m 


X{h,k) = - In 


{W{h,h)){W{h-l,l2-l)y 
{w{h - 1X2)) {W{WM-i)). ■ 


(151) 


If we assume an asymptotic mixed behavior of the Wilson loop of the kind {W{li,l 2 )) oc eyi\i[—KAhh — Kp2{li + k)], 
X converges to ka for large values of k and h. In our case, however, this asymptotic decay is meaningful only for 
h < Li/ 2 due to the periodic boundary conditions. In the ratio defining y, the perimeter contribution disappears 
because it is the same in the numerator and denominator, therefore the parameter y must go to zero in a deconfined 
phase (the ratio inside the logarithm tends to 1 and ka results 0), whereas it must be positive in the confined ones, 
corresponding to ka for large values of both k and k- 

As expected, the parameter y in the phases C and D converges fast to zero, thus showing the absence of an area 
contribution to the decay. 

For the phase D we considered, as an example, the point y = 5, z = 0.1. In this point xihik) < 5 x 10“^^ for 
Zi = 4 and Z 2 > 4 in a system size with Li = 8. The data show also a reduction by a factor of about 30 going from 
Zi = 2 to Zi = 3 and a similar reduction from 3 to 4. This is in precise agreement with the perimeter law shown in 
Fig. [T^and analogous data we obtained at smaller system sizes. 

The phase C presents a remarkable decay of the parameter y as well. For y = z = 1, the value of y drops to 
y « 5 X 10“^ for = 4 and Z 2 > 4 for both Li = 8 and Li = 9. Similarly to the previous case, y decays consistently 
when increasing the loop width. This behavior is therefore consistent with a perimeter decay and C behaves as a 
deconfined phase. 

Such a deconfined regime characterizing the phases C and D would be impossible in the compact QED due to the 
well-known results of [SIISHIIMI, and therefore we interpret the presence of these deconfined phases at t = 0 as the 
effect of the lattice structure of the PEPS which allows, for z —>■ 0, the decoupling of the 2D system into a collection 
of weakly coupled ID chains. 

The gapped phase B shows, instead, a totally different behavior. Here the Wilson loop is characterized by a faster 
decay dominated by an area law (see Fig. [12)3). In a system with width Li = 8, though, this area law is affected by 
strong finite size corrections for the loops with k = 6, 7 and k > 9, whose width is comparable with the width of 
the system. These finite size effects are totally absent instead for k < 6. The area law is indeed confirmed by the 
convergence of y to a value different from zero for Li = 8 and k < 6: Taking as an example the point at y = 1.32 and 
z = 1.77, y rapidly converges to the value y = 0.67 for all the values of li < 6 in (1511. We observe that in this phase 
the decay is much faster than in phases C and D. This leads to large numerical errors in the estimation of y for large 
loops {k > 10). For all the values with k < 6, 2 < ^2 < 10, though, y behaves like a constant, clearly indicating an 
area law contribution with > 0 (see Fig. 12 3). Phase B displays therefore a confinement of static charges. 

The data related to phase A are instead more difficult to interpret: For all the calculation at L 2 = 8,9 and 
2 < k < L 2 — 1, y presents an alternation in its sign dictated by the parity of the area; however, this cannot be 
considered a good order parameter for this phase, as the Wilson loops there obtain significantly low expectation values, 
which might lead to an accumulation of numerical errors in the calculation of y. Due to this reason, we checked also 
the value of the loops for smaller system sizes which reduce the decay rate: In this phase the results obtained for 
L 2 = 6 are very different from the ones of larger sizes and they present a positive limit of y without even/odd effect. 
If we consider, as an example, the point at y = 0.32 and z = 2.42 we obtain x{k = 2) « 1.5 and x{h = 3) « 0.08 for 

k > 10. 
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FIG. 13. Gap Z\ as a function of t for the given values of y and a and a cylinder circumference of L\ = 8. The gap does not 
close as fermions are introduced in the system, and according to the figure closes, presumably, closes for t —>■ oo. Inset: zoom 
of the data for small t 


From the numerical data it is therefore difficult to understand the nature of the phase A. We believe, however, this 
phase to be confining: if one expands the fiducial state of a single vertex around y = 0, 1/z = 0, an opposite behavior 
to that of the D phase is obtained: On top of the zeroth order vacuum, the leading O ( 2 ^”^) terms involve corners of 
flux loops meeting on a vertex. Thus the most significant nontrivial contributions in this state are excitations of single 
plaquettes (“glueballs” - the shortest possible excitations of the gauge field) which might imply that the A phase is, 
indeed, confining. 

A schematic plot of the phase diagram discussed above is given in Fig. El As a final remark on this topic, note 
that since we are dealing with the ground state with no static charges, all the possible states involve closed flux loops 
or infinite flux lines. 


2. The matter-gauge theory 


The introduction of a parameter t ^ 0 introduces the staggered matter fermions into the system as well. Most of 
our numerical results do not show any particular discontinuity in the gapped phases from the pure gauge case to the 
complete one (see, as an example, Fig. 13). However, it is still reasonable to expect a discontinuity in the t > 0 limit, 


as this point is where the fermions are coupled, or decoupled, from the gauge field. 

Our numerical data imply that the phases B and C, which used to be disconnected in the pure gauge model, 
seem to become adiabatically connected. The evaluation of the transfer matrix gap A still displays a minimum in 
correspondence of the critical line at t = 0, but, for t > 0, there is no evidence of a closing of this minimum for 
increasing system sizes. Therefore, the data seem to suggest that the critical line separating these phases in the 
upper panels of Fig. actually disappears as shown in the lower panels and it constitutes an isolated line in the 
three-parameter phase diagram, which may be an evidence for a discontinuity at t —>■ 0. 

Thanks to the introduction of matter, the set of non-trivial gauge invariant observables becomes richer and open 
bosonic strings delimited by fermionic operators provide further tools to examine the state {ipb)- In particular, due 
to the staggering of the matter fermions, we can distinguish two kinds of open string operators acting on the bulk of 
the system: 


hev 

=ipif n ^±(b)v^x. ■ 

bev 


(152) 


Afl represents the creation operator of a meson composed of a pair particle/antiparticle pair located on the even and 
odd lattice sites and Of and linked by a Wilson line defined on the path V. As in the previous definition of the 
Wilson loops, the signs must be chosen in order to fulfill the local gauge invariance: The Wilson line is oriented 
from the initial site to the final one Of in such a way that upward and rightward links are associated with 
whereas downward and leftward links with E_. J is instead a tunneling operator of a particle or antiparticle along the 
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FIG. 14. Expectation values of the meson operator (1) for various lengths I, with t — 1. Exponential decay is apparent in all 
the four phases. Data was sampled in the points A {y — 0, z = 5), B {y = 1.32, z = 1.77), C {y = 0.1, z = 0.1), D {y = 5, z = 0.1). 


path V, assisted by a Wilson line, and xy are both even for matter particles or both odd for antimatter particles. 
The signs +/— along the Wilson line must be chosen consistently with gauge invariance as well. 

We have numerically evaluated the expectation value of the meson creation operator {1) for a vertical meson of 
variable length I, with t = 1. The analysis was performed on a cylinder geometry with Li = 8 and L 2 = 40 + ^ with a 
starting point ei located at xn = 21. In all the gapped phases these expectation values decay exponentially with the 
length of the meson (see Fig. 14). This is consistent with a picture in which the transfer matrix Ts^., obtained by 
Eq. (141) with the substitution U = E+, has a non-vanishing gap between its two largest eigenvalues. 

Therefore the expectation value of alone is not enough to characterize the gapped phases, but it can be compared 
with the closed Wilson loops to provide information about the screening of the charges in our state. In particular, one 
can define the so-called horseshoe order parameter (closely related to the Fredenhagen-Marcu order parameter 
[721123]) associated with a rectangular loop on the lattice. Such an order parameter corresponds to the ratio of the 
(squared) expectation value of the meson operator on half of the rectangle C and the Wilson loop associated to 
C (see Fig. |T^ . Due to the limitations of the size of our system, we adopt a rectangle of size 4 x I (instead oi I x I in 
the ideal case presented in [7T1[72[71]). We define (the square root of) the horseshoe order parameter as 


P{1) = 




vWaO) 


(153) 


where the Wilson loop VF(4, 1) is associated to a rectangle C of dimension 4 x I with an odd I (in our staggered case), 
the matter particle and antiparticle associated to Ml are created in the middle of the two horizontal edges of C and 
the path V covers half of the rectangle C (see Fig. 15). Differently from the standard evaluation of the horseshoe 


order parameter, we are constrained to take one dimension of the rectangle fixed to /i =4 and thus we can only 
consider the thermodynamic limit in the vertical direction. 

p {1) is an order parameter for screening, assuming the Wilson loop has a perimeter law behavior, as expected for 
conventional gauge theories with a fundamental charge [ziii7a[7i|. Originally, it was formulated izaiTi] for the case 
in which the open edges of the horseshoe are in the temporal direction, and thus it is related to string breaking. For 
increasing dimensions of the considered loop, we can summarize the behavior of p in the following way: In a phase 
which is neither confining nor screening for the dynamical matter, p tends to zero; for confined or screened phases, 
instead, p tends to a finite limit. Indeed, in a deconfined phase free matter particles are strongly suppressed whereas 
the decay of the Wilson loop follows a slow perimeter law, therefore p falls exponentially to zero with increasing 
loop sizes 1. If, instead, the Wilson loop decay is dictated by a charge screening mechanism, the numerator and the 
denominator in (153) have a comparable decay and in the limit of large loops a finite value p ^ 0 will be reached 

[niizii- 

p may be associated to a line tension m between matter charges or even to the conductance properties of the 
matter [75]. Following dl, within a theory with a bosonic “frozen Higgs” matter field in the limit ^ 00 , p —>■ 0 

implies deconfinement while p 7 ^ 0 is a manifestation of a confining phase. In general, however, and in particular 
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FIG. 15. Schematic representation of the operators entering the horseshoe order parameter. The red rectangle C, with 
dimensions 4x3, represents the Wilson loop in the denominator of Eq. ( |153[ ). The dashed violet line depicts instead the 
meson creation operator in its numerator, associated with half of the rectangle labelled by V. We specify the bosonic operators 
involved along each edge of the loop, and the two fermionic operators (colored squares) entering the definition of MK 


for a theory with fermionic matter, such as in our case, p —>■ 0 provides only a tool to distinguish screening from 
non-screening regimes, while confinement itself cannot be deduced as its consequence. 

We have evaluated, for representative points at each region of the phase diagram for t = 1, both p{l) for various 
lengths {Li = 8, L 2 = 40 -I- 1) and the expectation values of Wilson loops with different width and length. The 
results may be seen in Figures [TT] andIt is indeed important to consider both order parameters: as explained, the 
horeshoe/Fredenhagen-Marcu arguments are valid when the Wilson loops follow a perimeter decay, which is expected 
in the presence of a fundamental dynamical charge in a conventional (e.g. Kogut-Susskind) theory. However, in 
our case, due to the truncation and the PEPS construction, one should not necessarily expect only a perimeter law 
behavior, and, indeed, our numerical calculations show a more complicated scenario. According to the numerical 
results, it seems that regions A, C, D obey a perimeter law for the Wilson loop, while the region B follows an area 
law (see Fig. 16). 

After having presented the Wilson loop results, we can go on to the results of the horseshoe parameter p. It may 
be clearly seen from the numerical data, that p ^ c 0 for increasing values of I, in the three regions A, C and 
D, whereas p —>■ 0 for the phase B. A perimeter law behavior of the Wilson loop is apparent in A,C,D, therefore 
one could most likely interpret these regions as phases which exhibit screening of the dynamical charge. This might 
suggest that the limit t —>■ 0 is not analytical indeed, as the inclusion of dynamical fermions may turn the pure-gauge 
deconfined phases C, D into screened ones in the dynamical matter case. 

Region B displays a less straightforward behavior. If it had a perimeter law for the Wilson loop, p —>■ 0 would mean 
that the charges in this region are not screened. However, since it manifests an area law behavior we cannot claim 
that with certainty. It might be still a confining/screening phase, due to the area law behaviour of the Wilson loop, 
although one must keep in mind that the area law is known as a valid confinement criterion only for static charges. 
Wilson loops with an area law in a conventional (Kogut-Susskind) theory with fundamental charges do not exist, and 
thus this result cannot be clearly interpreted. 

Another interesting feature of this phase diagram, as mentioned before, is the lack of a phase boundary between the 
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p=2(/i+y p=2{k+h) 

FIG. 16. Area/perimeter law behavior of the Wilson loop, for t = 1, at four representative points of the phase diagram. 



/ 



FIG. 17. Left panel: The horseshoe order parameter p for various lengths Z, with t = 1, for the four phases A, B, C, D described 
in the text. Data was sampled in the points A (j/ = 0, 2 = 5), B (y = 1.32, 2 = 1.77), C {y = 0.1, 2 = 0.1), D {y = 5, z = 0.1). 
Right panel: A schematic plot of the phase diagram for the gauge theory with dynamical fermionic matter {t = 1), with y,z > 0. 
The A, C, D phases seem to have a charge screening mechanism, while B’s behavior is unclear, but it might not screen. 


B and C phases, although they manifest different physical properties as described above. This may be understood, 
perhaps, as part of the discontinuity in the f —>■ 0 limit, where the phase boundary of i = 0 disappears. Such a 
crossover between the two regimes is also reminiscent of the non-separated phases in the Abelian models studied by 
Fradkin and Shenker m, although those models are very different from ours, since they involve bosonic Higgs field, 
rather than the fermionic matter discussed here. 


A schematic plot of the phase diagram is given in Fig. We further corroborated the structure of the phase 
diagram by calculating various expectation values (Wilson loop, ’t Hooft loop and meson string) for 2 = 1.5 as a 
function of y S [0,3]. The results are shown in Fig. 18 The expectation values display pronounced peaks and cusps at 
the points where the gap gets small, consistently with a possible non-analytical behavior in the thermodynamic limit. 
This indicates quantum phase transitions at the intersections of the line 2 = 1.5 with the phase boundaries shown in 
Fig. 17 On the other hand, a similar plot, of the expectation value of a 4x 3 Wilson loop along the line 2 = l/2 + 4/3y 
can be seen in figure [T^ This line goes, as y, 2 increase, from the C region to the B region, however, no non-analytical 
behavior related to a phase transition is seen, but rather two small, seemingly smooth peaks, presumably marking 
the edges of some transition region between the B,C regions. 
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FIG. 18. Magnitude of the expectation values of the Wilson loop (size 3x5), ’t Hooft loop (size 3x5), meson string of length 5 
and the gap A for = 6, t — 1, z = 1.5 as a function of y. The expectation values seem to behave non-analytically at around 
y ~ 0.65 and in the region 1.85 < y < 2.25 indicating transitions between the phases A and B and B and D, respectively. 
The Wilson loop and meson expectation values are plotted in a semilogarithmic scale (left axis), whereas the gap A and the ’t 
Hooft loop refer on the normal scale (right axis). 
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FIG. 19. The expectation value of a Wilson Loop along a cut through the regions B, C (the line 2 = l/2 + 4/3i/). A ’’transition 
regime” may be seen between the two peaks, however, they seem too smooth for a phase transition, agreeing with other 
numerical findings which refute the existence of a phase boundary between the regions. 


V. SUMMARY 


In this work, we have presented a method to extend the class of symmetric PEPS to local gauge symmetries. 
In particular, we have focused, as a first demonstration, on a 2 + 1 dimensional truncated compact QED with 
fermionic matter, but the methods presented throughout this paper are generalizable to other gauge groups [77], and, 
theoretically speaking, also to higher dimensions. 

We have shown how to construct, using the Gaussian formalism, globally invariant states for staggered fermions 
on a square lattice, satisfying fundamental physical symmetries; we have studied such states analytically, focusing on 
their parent BdG Hamiltonians and phase diagram. Then, we introduced suitable gauge degrees of freedom on the 
lattice links to gauge the symmetry and made it local. In this way we obtained a set of locally U{1) invariant states 
of both fermionic matter and gauge fields, with similar physical symmetries (and parametrization) as those of the 
purely fermionic, globally invariant states. 

Important properties of the PEPS have been discussed and demonstrated as tools for the study of the states we 
constructed. Eor example, we have used the spectrum of the transfer matrix of the PEPS, which is easily obtained 
numerically, to define the phase diagram of the system, by distinguishing gapped and gapless regions. We have also 
demonstrated the possibility of measuring several key physical observables for the study of the PEPS, such as Wilson 
loops or horseshoe order parameters. Through these results we have been able to sketch and characterize a phase 
diagram for truncated compact QED in 2 + 1 dimensions, which presumably includes both confined and deconfined 
phases in the pure gauge case and possibly some screening phases in the presence of dynamical matter. However, the 
main goal of this paper is to provide a demonstration and a proof of principle of the power of PEPS for lattice gauge 
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theory calculations, and the ability of PEPS to encode local gauge symmetries and describe the states of lattice gauge 
theories. 

The methods presented here can be further generalized and exploited for the study of other lattice gauge theories 
m With the application of efficient numerical techniques and methods, they may be used for a massive, systematic 
study of lattice gauge theories, and help shed light on their physics - e.g., such states, presumably with a larger 
bond dimension, may be used as variational ansatze for the Kogut-Susskind Hamiltonian [^. The road along this 
direction is still long, but with the use of current ideas and techniques, as given in this work and other ones recently 
published (see for example [351 [3B]), one could add tensor networks, and PEPS in particular, to the stack of available 
computational methods in high energy physics, and benefit a lot from the rich variety of options they suggest. 
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Appendix A: The Covariance Matrix Parametrization 

In this Appendix we investigate the properties of the Gaussian fiducial states \F) fulfilling the global U(l) gauge 
symmetry under the point of view of their covariance matrix. This analysis provides a deeper insight of statement 
and poses the basis for the explicit calculation of the parent Hamiltonian and correlations of the state \tjj{T)). 

A covariance matrix for (Dirac) fermionic operators aj , may be decomposed into the following block structure 

[78]: 


r = 


n Q 
Q 


where 


(Al) 


Qki = 2 


Tlu = - 


ak,al 


(A2) 


From these definitions, it follows that TZ is anti-Hermitian, TZ = —TZ\ and that Q is anti-symmetric, Q = — A 
pure state satisfies = |l. 

Our PEPS construction for the globally invariant state relies on the fiducial Gaussian states \F) defined by the 
Equations ( 20|21 ). Let us concentrate on an even vertex (the case of an odd vertex will be obtained by exchanging 
positive with negative modes). There, if we order the operators such that the negative modes come before the positive 
ones, we obtain the following block structure: 


F = 


( Faa 
TZba 

Qaa 

\ Qba 


TZab 

Qaa 

Qab 

TZbb 

Qba 

Qbb 

Qab 

TZaa 

TZab 

Qbb 

T^ba 

T^bb 

obtained 

from 


\ 


(A3) 


Let us consider T, the covariance matrix of the state obtained from the fiducial state by a gauge tranformation: 

(A4) 


From the definition of G (181 we obtain. 


F = 



|E) = 

\F). 

F-aa 

e^^^ab 


^'^Fba 

Fbb 

Qba 

^^*Qaa 

Qab 

Faa 

Qba 

e'^^^Qbb 



Qab_ ^ 

e-^^^bb 
e-^^^TZab 


(A5) 
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or simply 


r = v{4>)rv^ W 

where 


V (<P) = 


and 


Vo 


1 0 0 0 \ 
0-1 0 0 I 
0 0 -10 
0 0 0 1 / 


(A6) 


(A7) 


(A8) 


is the generator of the covariance matrix transformation correponding to the gauge transformation. 

Statement Al. F = F, if and only if the fiducial state is gauge invariant, i.e. G |J^) = q |F) for some integer q. 

Proof: 

1. G\F) = q\F) F = F: Since |F) is an eigenstate of the gauge transformation, |F) = \F) - the two 

states only differ by a global phase, thus the corresponding expectation values of all operators are identical, in 
particular the ones from which the covariance matrices are constructed. 

2. F = F G\F) = q\F): If the covariance matrices of both states are equal, and the related states are Gaussian, 

they correspond to the same state and thus may be described by equivalent density matrices (as Gaussian states 
are completely classified by their covariance matrices). Moreover, the states are pure, and thus may differ only 
by a global phase. Thus we deduce that If/) = \tp), or that |F) = \F). Since the spectrum of G 

contains only integer eigenvalues, we deduce that there exists an integer q such that |F) = \F), and 

that completes the proof. □ 

What, then, does the equality of the covariance matrices imply on the block structure? 

Statement A2. F = F, if and only ifTZab = 0,77ha = 0, Qaa = 0, Qbb = 0. 

Proof: F = F ii and only if [Vq, T] = 0 (considering an infinitesimal transformation). We denote 


Vo 


0 \ 

0 -Vo)^ 


% 


laa 0 \ 

0 -Ibb ) ’ 


where 1 labels suitable identity matrices, and we calculate the commutator: 


[Vo,r] = - 




(A9) 


(AlO) 


this will vanish as long as TZ commutes with Vq and Q anti-commutes with it, which ensures the block structure. □ 
Note that if \F) is gauge invariant, the expectation values of all the non gauge invariant operators with respect to it 
must vanish, implying the block structure: this implies one direction in the combination of the two above statements. 

Statement A3. If G |F) = q |J^), then Tr {TZbb) — Tr {TZaa) = —iq + | (Ap — A„). 

Proof: G\F) = q\F) means that 


{F\G\F) = q 


(All) 


Recall that 

{F\ G\F) = J2 {F\ bibk \F) - {F\ 4^^ \F) 

k k 


(A12) 



and 
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\F) 



ak,a 


t 

k 


1 |J^) = iTZkk + 


1 

2 


therefore, 


(F| G\F)=i (Tr - Tr (7^,,)) F^iNp- N„) 


(A13) 


(AM) 


which completes the proof. □ 
As a corollary, we may state: 


Statement A4. G\F) = q\F) if and only if the corresponding covariance matrix (AS) satisfies the following condi¬ 
tions: 

1. TZab = 0, TZba = 0, Qaa = 0, Qbb = 0 

2. Tr {TZbb) - Tr {TZaa) = -iq + f {Np - A„) 

Note that FF^ = jl as well, but this is a result of the purity of the state, independent of gauge invariance. The 
statement is almost equivalent to statement the only missing part is q = 0, which may be also proven using the 
covariance matrix approach (see below). 

The block structure is left invariant under canonical (unitary) transformations of the type TZ —> UTUA\ Q —> 
UQU^, as long as U is decomposed into 


U = 


Ur, 0 
0 Ur, 


(A15) 


i.e., U G U {Nn) X U (Np) and [U, Vq] = 0. These transformations are passive, and furthermore do not mix the 
positive and negative modes, which guarantees that the symmetry, and thus the structure of the covariance matrix, 
is conserved; it is trivial for 72., and straightforward for Q, as Q —>■ UQU^ results in 


Q 


0 

-UpQlM 


T 

P '^ab^n 


UnQabUp 

0 


Since \F) is a pure fermionic Gaussian state, it may be written in a BCS form |37j . 

1^) = n + Vkal'bl'^ |fl) 


(A16) 


(A17) 


where ul+vl = 1 (we make the assumption they are real, as shall later be clear from the singular value decomposition 


of statement A5), and the covariance matrices take the form (in the alternating ordering a, 6, a, 5,..) 


Qo = ®VkUkcry 0 0i 


(A18) 


72.0 


1 

2 


(1 - 2vl) 1 



(A19) 


where the uncoupled mode is present since we have an odd number of modes, 0i is the 1x1 null matrix, and li the 
the 1x1 identity matrix. 

The question then, is, whether this canonical BCS form is the result of a transformation U gU {Nn) x U (Np). We 
shall see that this is, indeed, the case: intuitively, this must hold for symmetry reasons. We now show this explicitly, 
using statement 


Statement A5. The Bogoliuhov transformation into the BCS form preserves the gauge symmetry. 


Proof: First, perform a singular value decomposition of T in Eq. (211, to obtain 


T = WrrKWl 


(A20) 
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where Wn is a unitary 7V„ x 7V„ matrix, Wp is a unitary JVp x Np matrix, and A is a Nn x Np diagonal matrix, whose 
diagonal (square) part includes the eigenvalues which are all real and non-negative. 

The canonical transformation U of the kind ( A15[), such that 


= (Hn) 


ij 


■ bi — (Wp)jj bj 


(A21) 


with 


Ur, = wl-, Up = w; 


brings the state into the (normalized) form 


\F)=Af i/^exp 


|f]) 


(A22) 


(A23) 


which is the desired BCS state. We identify 


Af ^ — Vf^. 


(A24) 


Using ul -I- = 1 we finally obtain 


Uk = 




Ar 


Vk = 






A? 


(A25) 


and 


AA = (1 + A 


(A26) 


From this, we can evaluate the BCS covariance matrices, TZq and Qq. Let us now write them in our usual ordering 
of the modes (unlike in equations (A18)-(A19)). In the explicit form below, we assume that iV„ = Np + 1 (even 
vertices), but the following arguments will also hold for Nn = Np — 1 (odd vertices). 


/ 


TZq = 


0 

2 l-l-A? 
0 
0 
0 

0 

0 


0 

0 

0 

i 

0 

0 

0 


0 

0 

0 

0 

1 1-A) 

2 l+A? 

0 

0 


0 
0 

0 
0 
0 

0 ^ 


0 

0 

0 

0 

0 

0 

2 I+AU 




(A27) 


Qo = 


/o 

0 

0 

0 

0 


0 

0 

0 

0 

zAp 


(TUA!) 

0 0 

0 0 


iAa 




0 

^Ai 

(l+A?) 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


(A28) 
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According to statement |A3[ Tr {TZhb) — Tr {TZaa) = —iq + \ {Np — Nn). From this, assuming that |-/Vp — 7V„| = 1, we 
immediately get that q = 0: the only states suitable for this representations are the ones with zero static charge, 
which agrees with statement □ 

If we denote the covariance matrix of the BCS modes by Iq, the covariance matrix F for the original modes contains 
the blocks 


'R = U^'RdA 

and 


(A29) 


Q=U^Q^U 


and can be obtained as 


r = 


0 

0 W 



(A30) 


(A31) 


As a corollary, we deduce that the most general gauge invariant fermionic Gaussian state with no static charges 
can be parametrized by either the set of real numbers and the canonical transformation U G U (Nn) x 

U (-/Vp), or the 5x4 complex matrix T. We shall adopt the second parametrization, and see how one may reduce the 
number of parameters when the desired symmetries are demanded. For that we shall now turn to the construction of 
the non-local state. 

The covariance matrix of the fiducial states, both even and odd, is the one previously investigated, up to the 
difference between even and odd sites (staggering). For even sites = 5,Np = 4, and the other way around for odd 
ones. We may deduce that both the even and odd sites possess the same covariance matrix, with a different ordering of 
the modes: in even sites, the negative modes will come first, thus having the ordering {ip, l+,r-,U-, d+, L, r+, u+, d-}, 
while on the odd sites, the positive modes come first, i.e. {ip,l-,r+,u+,d-,l+,r-,u-,d+\. 


Appendix B: The Gaussian Mapping 


Since we have a translationally invariant state, the PEPS may be calculated from the fiducial states and the bond 
states using a Gaussian mapping [55| . 

This is done in terms of Majorana fermions {cfc}: for any fermionic mode either positive or negative, one defines 
the Majorana operators 


C2fe-i = ak + a 


t 

k ’ 


C2fe = i 



(Bl) 


The covariance matrix of a state is proportional to the commutator of the corresponding Majorana operators: 


^Im = ^{[ci,Cm]) (B2) 

We denote the Majorana covariance matrix of the fiducial state of a single vertex, \F), by M. It is an 18 x 18 
real, anti-symmetric matrix, using the following convention for the mode ordering: { 1 p,l^,r-,l-,r^,u-,d^,U-^.,d-} 
for even vertices and {'0,r_|_, r_, (i_, u_, for odd ones. Note the positive-negative correspondence, in 

accordance with the staggering and the translational invariance. That means that the M matrices will take the same 
form on both the even and odd sublattices, but in two different bases. M will be used as the Gaussian channel, which, 
following [S3], may be decomposed into 

M = ( Jt o ) (B3) 


with A being the sub-block for the physical Majorana modes (a 2 x 2 matrix), D for the virtual ones (16 x 16) and 
B is the sub-block of physical with virtual modes (2 x 16). 

The bond states \H) ,\V) have the same covariance matrices: 

/ 0 ax 0 0 \ 

p _ -CTx 0 0 0 

° 0 0 0 

\ 0 0 -ax 0 / 


(B4) 



44 


where the ordering of the blocks is according to the ordering of the corresponding blocks in the covariance matrices 
M of the fiducial states, defined above. Due to the positive-negative correspondence, this means, for example, that 
the variance matrix of the first virtual mode with the second one is the identity of these “first” and “second” 
modes changes depending on the parity of the bond, whether it connects an even vertex to an odd one, or an odd to 
an even. However, the mathematical form is identical. 

The complete covariance matrix for the bonds will take the form (for an L x L lattice): 


= 00 


Perm (L, 1, 2,..., L — 1) i 


/ 0 CTj: 0 0 \ 

0 0 0 0 

0 0 0 (Ti 

\ 0 0 0 0 / 


+ Perm (2, 3,...,!/, 1) i 


/ 0 0 0 0 \ 

—(Jx 0 0 0 

0 0 0 0 

0 0 —Ox 0 y 


(B5) 


where Xi labels either the rows or the columns and Perm(ji, •■•iJl) is a permutation matrix whose nonzero elements 
are the entries The Fourier transform of Fin is simply: 


(^1 ? ^ 2 ) 


/ 

0 

_ 

(Jx^ 

0 

0 



/ 

0 


0 

0 

\ 


_ ^ — 

0 

0 

0 




0 

0 

0 


0 

0 

0 

_ 

(7xG 


0 


0 

0 

0 

_ 

(7x€ 


V 

0 

0 


0 

J 


1 

0 

0 


0 

/ 


(B6) 


We are now ready to apply M as a Gaussian channel on Fin. All the modes are aligned in a proper way, and 
since everything is translationally invariant we can work in the momentum space and obtain the momentum space 
covariance matrix of the state \'>l} {T)) [37]: 


“1 dT 


Gout (k) = A + H(D-Giu(k))-^H 


with A, B, D being the blocks in (B3). Following m, one obtains 

1 


Gout (k) = 


iPjk) Q (k) 
-Q(k) -iP{k) 


Vik) 


iPo (k) Qo (k) 

-Qo (k) -iPo (k) 


(B7) 


(B8) 


where 


I?(k) =det (D-Gin(k)). 


(B9) 


By writing the inverse in Eq. (B7) in terms of the adjugate matrix, it is easy to discern that the functions Po(k) and 
Qo(k) are trigonometric polynomials with a maximum order of 4 in ki and in ^2 individually (the same is true for the 
determinant I?(k)). However, since in the fiducial state \F) positive and negative modes are not mutually entangled 
due to the gauge symmetry condition, these polynomials might turn out to have lower maximum orders in practice. 

Furthermore, due to Goot(k) = —Gout(k), P (k) is a real function. Q (k) may be decomposed into its real and 
imaginary parts. 


Using Gin (-k) 
relations apply: 


Q (k) = i? (k) + il (k) 


(BIO) 


Gin (k), and similarly for Goutj as A,B,D are real, one obtains from (B8) that the following 


P (-k) = -P (k), R (-k) = i? (k), / (-k) = -/ (k) 

(Bll) 

and from the purity of the final state 


I?(k)=^i?2(k)+/2(k)+p2(k) 

(B12) 

Furthermore, since d (k) is real. 


X>(-k) =V{k) =X>(k) 

(B13) 

and thus 


Po (-k) = -Po (k), Ro (-k) = Ro (k), lo (-k) = -lo (k) 

(B14) 


hold as well. 
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Appendix C: Proofs of statements 5-7 


Some of the analytical results for the globally invariant case, namely statements 5-7, are strongly based on the 
PEPS construction, and thus their proofs have been omitted from the main text. We hereby give the full proofs of 
these three statements. But for that, we begin, first, with a formulation of the PEPS in momentum space. 

Eollowing the convention for physical fermions (79), we define a similar Fourier transform for virtual fermions: 


at = 


bi = 


1 


'jLiL-, 






y/LiL- 




(Cl) 


k- 


To manifest translational invariance better, we swap the a and b operators on odd sites, i.e. exchange all the a 
operators there by b operators and vice versa. This results in a uniform expression for the operator A everywhere, 

A (x) = exp (C2) 

where the j indices are 1, ...,4, involving only virtual fermions, while these of i are 0, ...,4, with qq = ip the physical 
fermion. 

The state \ip{T)) defined in ([^ involves a product of these operators everywhere, which results in a summation, in 
the exponent, over all the sites; we can perform a Fourier transform, and obtain 


(x) = (k') A (0,0) A (tt, tt) A (tt, 0) A (0, tt) (C3) 

X k' 

where k' involves only half of the k values in the Brillouin zone (the “positive momenta”), except for (0,0), (tt, tt) , (tt, 0), (0, tt), 
and, as a result of the Fourier transform, 

A (k) = exp (Tijuj (k) b] (-k)) exp (t.^oI (-k) b] (k)) (C4) 

We continue with Fourier transforming the projectors. The product of all projectors may be denoted as 

\{HV}){{HV}\=(^\H){H\®\V){V\ (C5) 


(note that as these operators are quadratic in the fermionic operators, a tensor product is well defined) with 


\{HV}) = n 0) ^) K 0) 

k' 

and 

k)) X exp {Sij (k) b, (k) bj (-k)), 


0 0 \ 

0 0 0 
0 0 ■ 

0 0 

Since (fl„| {HV}) is just an irrelevant constant (equal to 1), the physical state is simply 

\iij{t))= 

and we can decompose it into a tensor product of momentum states, 

li’ (T)) = (^ IV' (k')) C) \ip (0, 0)) \ijj (tt, tt)) \ip (tt, 0)) \ip (0, tt)) 


B (k) = exp {Sij (k) (k) (- 


where 


( 0 


s (kx, ky) = 


iki 


0 

V 0 


(C6) 


(C7) 


(C8) 


(C9) 


(CIO) 
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where we separated paired and unpaired modes. The paired ones have the form: 

|V^(k')) = (f^„|S(k')4l(k')|r!) 


(Cll) 


with | 0 ^), |r 2 ) being the vacua of the relevant momentum subspace (virtual and general respectively), and a tensor 
product structure is again well defined as all the states in the product have an even fermionic parity. 

For k = (0, 0), (tt, tt) , (tt, 0), (0, tt), one defines instead S (k) = S (k) /2: the operators B are defined using the S 

matrices instead of the S ones, and A (k) = exp (Tija] (k) (k) j . 

Proof of statement^ Denote the blocks of the covariance matrix of the fiducial state by A, B, D, the initial bond 
state by \ipin) and the output state by |'0out)- Then (assuming the number of virtual fermionic modes is even, which 
is, indeed, our case), the norm of the output state is given by [53] 


(V'outIV’out) = Pf - Gin^ (■i/'inIV’in) 


(C12) 


where Pf denotes the Pfaffian. 

In our case, the fiducial state was not normalized, and thus the fact that |' 0 out) is not normalized is not a mere 
outcome of the Gaussian mapping. Therefore we modify this formula by multiplying by a further constant 7 . That 
shall be done carefully, considering the paired and unpaired momenta separately. For the paired momenta, D = D(BD 
and Gin = Gin (k) 0 Gin (—k), while for the unpaired ones, D = D and Gin = Gin (k). Next, since the fiducial states 
lack some normalization factors, we deduce that such factors should appear twice for paired momenta, and once for 
unpaired momenta, i.e.: 


{if (k) 1^ (k)) = yPf [D 0 D - Gin (k) 0 Gin (-k)] (V-m (k) |^i„ (k)) 


for paired momenta, and 


(V- (k) IV' (k)) = V^Pf [D - Gin (k)] (V'in (k) IV-in (k)) 


for unpaired momenta. 

For the paired momenta, since V (k) = T) (—k) = det {D — Gin (k)), one simply obtains: 


V{k) = 


{if (k) |V:(k)) 

l{lpin (k)l^in (k)) 


whereas for the unpaired ones: 


V{k) = 


|a(k)|V|/3(k)|^ 

2567 


i5(k)r 


{if (k) Ilf (k)} ^_ 

7 (V'in (k) IV'in (k))^ 2567 


(C13) 

(CM) 

(C15) 

(C16) 


thus we obtain continuity if a (k) = a (k) - which, as we shall show, holds. Given this assumption which will shortly 
be proven, we define 

E{k) = \a{k)f + \P{kf (C17) 

as the dispersion relation we work with. This proves statement □ 

We shall now turn to the calculation of the functions a (k) ,j5 (k) - the proofs of statements and For that, we 
first prove a more general statement. 

Statement A6. The expectation value 

F {A, B,C, D) = (f 2 | (k)h^(-k)gD^a(k)a((-k)f)t(k) (C18) 

is given by: 

F (A, B, G, D) = det (ADBGT 0 1) (C19) 

Proof: Let us calculate F explicitly. The exponentials have to be expanded, but since the creation and annihilation 
operators must be balanced, they should all be expanded to the same order, and thus one obtains 




Ai^j^.-.AiNjKBkih X 


G„ 


N 

,...G, 


(C20) 


7 ~) 7 ~) yii...iNji...jNki...kNh---lN 

OtN0N-^7lSl'---^-fNSN ■^Oci...OtN01.../3N7l---yNSi...SN 
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where Z is the vacuum expectation value of the product of creation and annihilation operators, which we shall now 
calculate: 




W Oji (-k) (k) (-k) X 

fefci (k) hi^ (-k) (k) bi^ (-k) X al^ (k) (-k) (k) b^^ (-k) x 

«7l (^) (-J^) C 1^) ■ 


(C21) 


After changing the operators’ positions, in a way which leaves the sign invariant, we can decompose this term into 
a product of four vacuum expectation values. 


(k)(k) (k) (k)) x 

(oji (-k) (-k) a\^ (-k) (-k)) x (bk, (k) ...6fc« (k) (k) (k)^ x 


bi, (-k) ..Mr. (-k) bl^ (-k) ...b\^ (-k)) = (5^^^ 


*iv xji...jrr ski...krf di...l]v _ 

Q;iv...ai^7iv--.7i^5jv---i5i ^/ 3 n .../ 3 i 


(C22) 


Hi...in xji...jN rfei-.-few 
^ai...ajv^7l---7-N^i5i...i5N ^Pi.../3n 


where the Kronecker symbol is defined by = ±1 if *i, are a cyclic/anti-cyclic permutation of ai...aAr, and 

0 otherwise [THUD], and the last equality results from exchanging the order of the lower indices in all the Kronecker 
symbols, which leads to the same sign in all four of them and thus leaves the total sign invariant. 

Note that i® antisymmetric under an exchange of any two j indices: let us demonstrate that 

by exchanging ji,j 2 - Then, 


A. . A. . A. . — A . - A. . A. . A. . A. . 

'^ai...aN “ '^ai.. .a r; '' ’■JVj ]V 

An important property of the Kronecker symbol is that if is antisymmetric m then: 

JLxu...iNA^. . = M- 


(C23) 


(C24) 


and thus, using the antisymmetric behavior of ' 7 , ...Ai^j,.j and similar symbols, we can significantly simplify 

F (ICM: 


^ = E 


N 

Cn,, ft, ...C, 


1 

M 


Hi. ..In 4. . A n,, R, , V 

^ai...aN^~ii...^N^&i...&N ^/3i‘ ■■■k^kr.lr. ^ 


ai;fli ■■•'- 2 q:jv/ 3 jv-^ 7 i^i ■■•-^ 7 Af' 5 iV ~ ^ ^ 

N 


1 

m 


Hi...in s:ki...kN Hi...In ^ 


{-^^)iiSi ■■■ ^kili---Bkr.lr.C!aiPi---Car./3N — 

AT \ •/ 

E (... . 


(C25) 


N 


One may also use the Kronecker symbol for the calculation of the characteristic polynomial of a matrix M |80j . 

(C26) 


PM {x) =det {M-xl) = Y^ i-x)^ 

N ^ ■/ 


and thus we finally obtain that 

F (A, B, C, D) = padbci (-1) = det {ADBC^ + 1) = 

= det {DBC'^A + 1) = det {BC^AD + 1) = det {C^ADB + t) (C27) 


where the many possible forms are possible thanks to the symmetry properties of F. □ 
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Now we can finally turn to the proof of statement^ We start with the calculation of a. Note that 

a (k) = (flpIV') = (f^l B (k) A (k) \n) 

where jn) is the total, both physical and virtual, vacuum of the relevant momentum subspace. 

Using the definition of the state |'0 (k)), one can simply write a as 


(C28) 


(C29) 


Since there are no physical annihilation operators, the physical part of A does not contribute, and one may replace 
the vacuum by the virtual vacuum, and T by r. In this form, a (k) = F (5", 5", r, r). Then one simply obtains 


a (k) = det {St'^ St + 1) 


(C30) 


whose straightforward calculation leads to Eq. ( |89[ ). 

Similarly, 

/ 3 (k) = ( 0 | l/; (-k) V'(k) e‘S'o«dk)a,(-k)gSMfcfc(k)6i(-k)gT„;5at (k)fet (_k)gT.,sa+(-k)ht (k) ^ 

— ('.Ol p‘^^’j“hk)aj(-k)pSfeif)fc(k) 6 i(-k)pT„^at (k) 6 +(-k)pr^.5ot (-k)bt(k) _ 


lim ^ (U| e' „ „ 

X —^oo .X. 

lim ^FiSx,S,T,T) = lim 4 det + 1 ) 

X^aoX X^aoX 


(C31) 


where Sx = —© S (li is the 1x1 identity matrix). This results in ( [90| ) and completes the proof of statement 

□ 

The only remaining task, which shall result in the proof of statement [Tj is the calculation of a for the unpaired 
modes, k = ( 0 , 0 ), (tt, tt) , (tt, 0 ), ( 0 , tt). 

5(0,0) = |n) (C32) 

with all the operators at k = 0, and similarly for a (tt, tt). We expand the exponentials such that all the creation and 
annihilation operators are balanced, and obtain 


^(0,0)=^ 


N 


1 


1 


Nj {2N)\ 




(C33) 


oq Qj, .. .a,„ bk, hi^-hk^bi^a^^ b\^. 


U hX 


Note that S'( 7 r, 7 r) = —S'(0,0), but as the S matrices always appear an even number of times in the expansion, 
5 (tt, tt) = 5 (0,0). 

Here, only iV = 0,1, 2 contribute - higher orders vanish due to the fermionic statistics. The zeroth order is trivially 
1. The first one is 


5i (0,0) — 2 bia^Sja^) ((i;/3i^fc/32 ^//32 )^ct2/52 — 


- 2Tr 

And the second, last order is 


(sWStt) =-2 {y^ + z'^). 


1 1 


0^2 (0,0) — ^ ^|l5Q\‘(i2^a|a4'^/3(/32/l3^4‘^»iJi'^*2i2'^fclh‘^fc2b'^ai/3iT'Q2/32'©3/33'©«4/34- 


(C34) 


(C35) 


Since the number of possible indices is 4, the Kronecker symbols may be decomposed as = e*^'^^*^'^ 2 ea 4 a 2 a 30!4 


(e. 


ia 2 (y. 2 ,(y. 4 ^ 

Q. . c. . 

Ijl'i2j2 ^'i‘2j2 


etc., 9-Ild. tll 6 n — dot (t) 3,S W 6 ll 3-S 

second order is simply det (r) = (jj"^ — z^Y. Altogether we obtain 

5 (0,0) = 5 (tt, tt) = 1 - 2 (y^ + zY + {y^ - zY"^ = [l - (y + zf'^ (l - {y - z)^) . 
A similar calculation for the other unpaired momenta yields 

5(7r,0) = 5(0,7r) = 1 - 2 (y2 - z^ + {y'^ - z^Y = (l - [y'^ - z^Y ■ 


= 1 , hence the 


(C36) 


(C37) 


Note that the right limit is obtained for the unpaired momenta k: indeed, a —> 5^ there, as a straightforward 
consequence of Eqs. (891,(90). This proves statement □ 
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